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ABSTRACT 


An autonomous spaceborne gravity ^radiometer mission is being 
considered as a post Geopotenlial Research Mission project. The 
introduction of satellite gradiometry data to geodesy is expected to 
improve our solid earth gravity models. This study explores the 
possibility of utilizing gradiometer data for the determination of 
pertinent gravimetric quantities on a local basiB. The analytical 
technique of least squares collocation is investigated for its usefulness 
in local solutions of this type. It is assumed* in the error analysis, that 
the vertical gravity gradient component of the gradient tensor is used 
as the raw data signal from which the corresponding reference 
gradients as'e removed to create the centered observations required in 
the collocation solution. The reference gradients are computed from a 
high degree and order geopotential model. The solution can be made in 
terms of mean or point gravity anomalies, height anomalies, or other 
useful gravimetric quantities depending on the choice of covariance 
types. Selected for this study were 3O'x30' mean gravity and height 
anomalies. Existing software and new software are utilized to implement 
the collocation technique. It was determined that satellite gradiometry 
data at an altitude of 200 km can be used sucessfully for the 
determination of 30'x30' mean gravity anomalies to an accuracy of 9.2 
mgal from this algorithm. It is shown that the resulting accuracy 
estimates are sensitive to gravity model coefficient uncertainties, data 
reduction assumptions and satellite mission parameters. 
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I. INTRODUCTION 


Gravimetric geodesy is in the midst of a period of rapid growth. 
During the last twenty~fivw years, geodetic information has grown 
formidably, primarily due to the application of space technology. With 
proper guidance and funding, this trend will continue at a pace 
unimaginable to those scientists involved at the beginning of the space 
age. Currently, many countries have offerred support for the rapid 
growth of geodetic science. In the United States, geodetic progress is 
made by many government', agencies, universities and private 
corporations who are involved in a wide variety of projects. Perhaps 
the largest civilian endeavor is one under the auspices of the National 
Aeronautics and Space Administration known as the The Geodynamics 
Project. The ob^ctivos of this project, as outlined in NASA [1983], are: 

"(1) to contribute to the understanding of the solid Earth, in 
particular the processes that result in movement and 
deformation of tectonic plates; and 

(2) to improve measurements of the Earth's rotational dynamics 
and its gravity and magnetic fields", 

These two objectives are inextricably related and perhaps a third 
objective could be included; 

(3) to form a useful geophysical information data base for the ubo 
of studies in advanced Earth geophysics as well as absolute 
and comparative planetology. 

Thus the primary objectives could be reformed into the single sentence: 
The primary goal of The Geodynamics Project is to formulate solid-Earth 
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models which more accurately depict the geophysical nature of our 
planet. 

Several types of observable phenomena exist which contribute to the 
formulation of solid-Earth models. Two such phenomena are the 
geopotential fields, the gravity and magnetic fields of the Earth, In 
this study, the gravity field is highlighted as the phenomenon of basic 
interest. An interesting number of by-products, or applications, con be 
made from the study of the Earth's gravitational potential (or, for short, 
the geopotential), These applications are nicely summarized in the NFC 
publication, Applications of a Dedicated Gravitational Satellite Mission 
(National Academy of Science [1979]). 

The rationale for improving and increasing the knowledge of the 
Earth's geopotential comes from many sources. Some examples include: 

Geophysics: Gravity information is used to infer isostatic information to 

supplement and constrain tectonic models. Geoid 

undulations can be correlated with density anomalies to 
constrain mantle convection models. 

Oceanography; Highly accurate geoid models when used with altimetric 
data yield precise knowledge on the sea surface 
topography which is useful in the determination of local 
ocean circulation patterns. With altimeter measurements 
over a period of time, a fourth or dynamic dimension can 
be included which may improve the understanding of the 
meteorological/ocean circulation interface. 

As can be seen from these few examples, improvements in gravity field 
knowledge have strong repurcussions within the geophysical community 
which can ultimately contribute to a unified earth model. There are 
certain geodetic applications which would benefit from improved 
gravitational potential knowledge. The major impetus stems from the 
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need to improve the accuracy of artificial satellite ophemerides. Many 
NASA projects being considered rely heavily upon how accurately a 
satellite orbit can be determined. Altimetry projects (such as TOPBX), 
mapping projects (such as MAPSAT), and other remote sensing missions 
are some examples of missions requiring precise orbits. 

With this brief introduction* some examples supporting the rationale 
for gravity model improvements have been mentioned. It is then 
worthwhile to describe a few of the projects currently under 
consideration to directly improve our gravity field knowledge. 

NASA is currently considering two types of instruments to be used 
in a gravity field measurment campaign; Satellite-to-Satellite tracking 
(SST) techniques are to be used in the Geopotential Research Mission 
(GRM); and satellite gravity gradiometry utilizing cryogenic technology is 
to be used in a follow on mission. The Satellite-to-Satellite tracking 
instrumentation is described in APL [1983] and a proposal for recovering 
geopotential information is discussed in Colombo [1984]. Satellite borne 
gravity gradiometers (the instrument of primary interest in this work) 
are summarized in Wells [1904]. Instrument descriptions and 
gradiometry theory will be discussed in the next chapter. 

It should be stressed that the two types of instrumentation (SST 
and gradiometry) are not competitive with one another. Rather, they 
are complementary since they are individually sensitive to different 
wavelengths of the gravity field, This is evident in Figures la and lb 
which were generated by the rapid error analysis procedure described 
in Jekeli and Rapp [1980], The graph illustrates mean gravity anomaly 
accuracy estimates from the six me .th GRM mission as well as those from 
a. six month gradiometry missions. It can be seen in Figure la that 
higher resolution is achieved when the sensitivity of the instrument is 
increased. Higher resolution is also achieved when the satellite altitude 
is lowered (Fig. lb) such as proposed in the space shuttle tether system 
(NASA [1984]) which would allow for altit d ' near 140 kilometers. Thus 











if the GRM project is fully supported, then the follow -on gravity 
gradiometry miaaion, if it ia also supported, will be under strict 
instrumentation and altitude requirements in order to make significant 
improvements in spatial resolution. 

In this work Ve prospect of utilizing local (as opposed to global) 
gradiometry data ,>v the estimation of mean gravity anomalies and 
height anomalies is explored. The theoretical foundation is supplied by 
least-squares collocation theory. The idea of using only local 
information was advanced due to the substantial computer effort 
required in global solutions. The purpose of this work is to test a local 
collocation algorithm for its accuracy, dependability, and ease of 
operation. 

The next chapter introduces the satellite gravity gradiometer (SGG) 
concept along with the units and coordinate systems involved. Mentioned 
too, are the design concepts, likely mission profiles, and a review of 
previous reduction methods. In chapter three, the general theory of 
collocation as applied to the local gradiometry reduction is expounded. 
The use of a high degree and order geopotential model will be explained. 
Chapter four presents the results of the error analysis and ensuing 
discussion. The final chapter attempts to provide a perspective of this 
work by drawing conclusions and suggesting alternative lines of 
investigation. 


2. REVIEW OP GRADIOMETRY, ITS INSTRUMENTS AND METHODS 


2.1 Concepts and Theory 

The baaic theory behind gradiometers haa been known since the 
early daya of torsion balances, when around 1900 a Hungarian physicist, 
Roland Etttvtis, developed a working variometer. Torsion balances had 
wide spread use during the 1920*8 for geophysical exploration. Their 
use declined due to the introduction of the more economical gravimeter. 
However, torsion balance measurements are still being made today but 
they are rapidly becoming obsolete with the influx of high technology 
gradiometers and inertial systems. 

♦ 

In order to intuitively appreciate the concepts underlying 
gradiometry, it may be worthwhile to discuss the theory, units and 
phyaical interpretation behind this system. A gravity gradiometer is 
capable of sensing and recording the values of the spatial second 
derivatives of the gravitational potential. The discussion begins by 
recalling that the potential ie a scalar quantity having Si-units of 
m a s -a . The first spatial derivative of the potential is generally 
associated with gravitation which is an acceleration quantity or vector 
quantity having units of ms“ a . Acceleration can be intuitively described 
as the rate of change of an object's velocity. The second spatial 
derivative of the potential is also a vector quantity which has units of 
s~ a . These units may seem physically meaningless at first, but it is 
easy to recognize that the second spatial derivative allows for a 
description of the spatial variation of the gravitational acceleration. For 
instance, the unit most often associated with these second derivatives is 
the Etttvtts (denoted E). One EtttvttB is 10“* s" a which can be thought 
of as 10"* ms~ a /m. In words, this means that one Edtvtts is a change of 
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10“* ms"* over a one meter distance. In more familiar gravity unite, 
one Etttvtts ia equivalent to a 0.1 mgal change over one kilometer. 


The nine epatial aecond derivativea of the gravitational potential 
conatitute the "gravity gradient tenaor" which ia conveniently expressed 


in matrix form by 



V,, 

v, a 

Vis 


v a , 

Vji 

v a . 


V,i 

V aa 



(II. 1) 


where V denotea the gravitational potential of a body and the aubacripta 
refer to the apatial derivativea (eg. V tI = »*V/*xi, Vj a s **V/*x,*x a , 
etc.). The 1, 2, 3-subscripts refer to directions aaaociated with an 
orthogonal coordinate system. As long as the directions remain 
orthogonal, then the gravity gradients along these directions must 
satisfy Laplace's equation, 

I Vu = 0 (II. 2) 

i=i 


that is, the diagonal elements of the gravity gradient tensor are subject 
to Laplace's condition. 

The formulation of the potential by spherical harmonicB begins by 
defining a geocentric coordinate system. Let the X 3 -axis of a geocentric 
Cartesian coordinate system coincide with the mean rotation axis of the 
earth. Next, let the Xt-axis lie in the mean equatorial plane as well as 
in the plane formed by the Greenwich meridian. The X a -axis completes 
the right handed Barth fixed Cartesian system by lying in the mean 
equatorial plane in the direction of longitude rr/2. The geocentric 
spherical coordinates (♦, X, r) follow immediately which are defined as 
shown in Figure 2. The earth's gravitational potential can then be 
written in spherical harmonic form as 



Figure 2. Geocentric cartesian and spherical coordinate systems 


V(», X, r) = ® [i + } ( S J* YU (CimCosaX + S| UB sinBX)P nB (sin«)] 

n=a »=o 

(II. 3) 

where C M , S n « are the normalised potential coefficients, Pw(sini) are 
the normalized associated Legendre functions, a is a mean equatorial 
radius value for the earth (e.g. 6378137 m for the 1980 Geodetic 
Reference System) and r is the radial distance of the computation point 

f at latitude ♦ and longitude X. Similarly, an expression for the reference 

$ 

| gravitational potential can be written as 
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U(4, A, r) = ® [i + YZ ( f )" E (C^osnA + Sn a sinaA)P na (>int)] 

n=* a=o 

(II.4) 

where the priaes denote reference potential coefficients which are 
approx iaat ions to the C ^ . S M in equation (II. 3). 

The disturbing potential, in the present notation is, in the well 
known formula, 

T (♦, A, r) = V(f, A, r) - U(t, A, r) (II. 5) 

The disturbing potential can also be written in a spherical harmonic 
expansion as: 


£ 

I 



$ 

\ 

V 


T(4, >, r > = “ [ E ( ? )" E ( eCjjaCosaA + eSnasinaAjPn^sin#) 

n=a a=o 

■ n 

+ ) ( f )" ^ (CimcosaA + S nB sinaA)P nB (si,nf)] 

n=a B-o 

(II. 6) 


where here the eC nm and eS^ terns are the discrepancies between the 
coefficients of the true field and the coefficients of the reference field. 
That is, more explicitly, 



(II. 7) 
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These discrepancies are better termed "coefficient errors” which implies 
that the reference (or "model”) coefficients are in error with respect to 
the unknown "true” coefficients (Colombo, [1980]). Since the coefficient 
errors cannot be found directly, they are approximated by the a 
posteriori variances for the computed model coefficients. These irrors 
will play a significant role in the covariance computations discussed in 
Section III.4. 

Next, the spatial derivatives of the various potential expressions will 
be examined. A connection between the well known principles of 
gravimetric geodesy and the spatial derivatives will be highlighted to 
provide an intuitive framework. 

The first spatial derivatives of (II.3) can be related to the 
components of gravity by (Heiskanen and Moritz [1967]) 


«x, 


= V 


x, + 


r x. * 


tx a 


+ f 


x a . 




= Vv + *y 




(II. 8) 


where 4 denotes the centrifugal potential. Similarly, the components of 
the reference (or normal) gravity are; 


Tx, = Ux, + *x,* 7x a = Ux a + *x a » 


7x, = u Xo + *x, • 


(II. 9) 


The gravity disturbance, 3 is defined as the difference between the 
true (neasured) gravity, | and the reference (node led) gravity 
(Heiskanen and Moritz [1967], section 6-1) 




(II. 10) 


Then, in geocentric Cartesian coordinates, the gravity disturbance 
components are (since the centifugal potential terms cancel): 



•V 

#u 

r ixi~ 


*T 

Xk 

" *xr 

•X, 

7x» 

‘ *x 


•V 

«u 

= <Xa~ 


•T 

Xa 

■ 'X, 

«x a 

Txa 

•Xj 


♦V 

«u 

* «X3- 


•T 

xs 

“ 'X," 

*X S 

7xs 

" *X: 


(II. 11) 


It is possible to differentiate each equation in (11.11) with respect to the 
three directions thus yielding the nine anomalous tensor components 
referred to the geocentric cartesian system 


Tm 


t;, = 1 1.^.) 

T,a = TXj^x,) 


t;, 

= 7x, (fi ’ c a ) 

laa ~ *X a ^ x a ) 

Tas = '?X a (tfx a ) 

(11.12) 

Tji 


Tsa = ‘5x a (fix 3 ) 

Taa = 7x, (tfx , ) 



where the primes denote that the gradients are referred to the 
geocentric coordinate system. This equation illustrates how the gravity 
disturbances are related to the anomalous gravity gradients referred to 
the geocentric cartesian system. The anomalous gravity gradients are 
then the spatial rates of change of the gravity disturbances. 

It is expected that the gravity gradient observation will be referred 
to a local level cartesian coordinate system. In this coordinate system, 
the x 3 -axis coincides with the radius vector in the spherical case (Pig. 
3). The x 2 -axis is oriented in a northerly direction and the X|-axis in 


% '% 
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2.2 Instrumentation 

For the sak ' of comprehensiveness, a brief overview of satellite 
gravity gradiometer inetrumentation is provided in this work. The 
discussion is by no means complete but only serves to increase the 
reader's awareness of the measurmont principles involved and the 
accuracies obtainable by the proposed instruments. References will be 
made throughout to aid the reader's desire for details. 

At present, two very good reviews of instrument theory and design 
exist. The first one by Forward [1974], presents a good historical 
perspective of gradiometer mission and system development up to the 
date of publication. Several types of designs are presented ranging 
from quartz balance systems, vibrating string systems, to rotating 
resonant torsional gradiometers. Since the review by Forward, 
applications of new technologies have shown remarkable improvements in 
gradiometry design. These new generation gradiometer systems are 
reviewed in Wells [1984]. This section attempts to highlight some of the 
designs mentioned in this more recent publication. 

In the United States, there are a number of research groups 
investigating gradiometer instrumentation (Table 1). Some of the groups 
have experienced continuous funding from various agencies for the 
purpose of development, others have had their funding cut-off. Thus 
this overview will contain only those groups known by the author to be 
actively funded. 

Basically there are two types of gradiometry sensors; Conventional 
and Cryogenic. Conventional sensors operate at room temperature and 
may have a higher level of maturity than cryogenic sensors which 
operate at very low temperatures (less than 4.2* Kelvin, the boiling 
point of liquid Helium at one atmosphere). Only two of the organizations 
listed in Table 1 are developing conventional gradiometers, these are 
Bell Aerospace and Hughes. 


The Bell Aeroepace gradiometer utilises a slowly rotating platform on 
which are mounted four accelerometers (Pig. 4). The gravity gradients 
are measured by differencing the outputs of opposing accelerometers. 
Due to the rotation of the platform, the gradient signal is sinusoidally 
modulated with a period at twice the rotational period. The use of four 
accelerometers facilitates the detection of minute gravity gradients by 
utilising "common mode rejection of the large acceleration" (Metzger 4 
Jincitano, [1981]). An orbital mission system has been proposed by Bell 
which utilises the same instrument concept. Miniature electrostatic 

accelerometers (MESA’s) are to replace the current operational 

accelerometers (the Mark VII). It has been reported that the MESA 

based gradiometer will have much lower noise levels than the current 
operational system (the MESA based system has an estimated noir.e level 
of around 0.035E (Wells [1984], p. 32)). 

The Hughes gradiometer is also a rotating instrument but the system 
design is quite different from its Bell gradiometer. The observable 

quantity ia the strain due to the torsional flexure experienced by the 
resonant cruciform mass-spring system (Fig. 5). The torsional flexure is 
coupled to the differential torque experienced between the two arms. 
Since the system is rotating, the differential torques, due to the 
gradients of the gravitational field, ere excited at a frequency twice 
that of the system rotation frequency. The differential torque, AT, is 
related to the gravity gradients through the following expression 


AT = [(Vxx ~ Vyy)cos 


2wt + 2V xy sin2<t] 


(11.14) 


where all terms are defined as illustrated in Figure 5 and where the x 
and y directions lie in the plane of rotation. The system noise level 
goal is 0.01E (la) using a 35 second integration time. It was not 


Table 1 

Organization! Active in Gravity Gradiometer R k D. 


Organization 
Bell Aerospace 
Bendix/Stanford Univ. 

Hughes 

SAO/PSN 

Sperry Defense Syst . 
Univ. of Maryland 


Gradioaeter Type 

Rotating Acceleroaeter 
Superconducting Cavity 
Oscillators 
Rotating Gravity 
Gradioaeter 
Gravity Radiation 
Cryogenic Levitated 
Balance Arn 

Superconducting Accelerations 
with SQUID readout 


Principle 

Investigators 

Metzger/ J inci tano 

Reinhardt/ 

Turneaure 

Forward 

Hastings 

Paik 


spin 

axis 



spin 

axis 


Figure 4. 




Operational concept of the Bell Aerospace gradiometer 

(from Wells, 1984, p.31). 
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reported In Welle (1984) whether or not the system noise goal has been 
achieved. 

The remaining or gar '.nations listed in Table 1 are developing 
cryogenic based gravity gradiometars, these are; Bendix/Stanford Univ., 
SAO/PSN, Sperry Defense Systems, and Univ. of Maryland. Brief 
descriptions and instrument sensitivity goals make up the remainder of 
this section. 

The Bendix/Stanford University gradiometer utilizes superconducting 
Cavity Oscillator (SCO) technology developed at Stanford. The SCO is 
capable of detecting minuto displacements on the order of 10"* 9 cm 
which it then converts to frequency form. Thus the output signal 
would be a frequency shift corresponding to the displacements 
experienced by a mass-spring accelerometer which are due to the 
gravity gradients. Bendix Field Engineering is developing the design of 
the SCO Based gradiometer. Their gradiometer, called the canonical 
gravity gradiometer, has six 3-axis SCO based accelerometers placed at 
distances of 4/2 from the origin of a cartesian coordinate system (Figure 
6). All nin components of the gravity gradient tensor can be 
approximated directly from the outputs of the 18 component 
measurements from the accelerometers. Recalling the definition of an 
Eotvos, the gravity gradient can be determined simply by differencing 
two accelerometer outputs and dividing by their along axis distances 
(4/2 or 4 in the canonical case). The noise sources affecting the SCO 
based gradiometer have been investigated and it has been shown that a 
gradiometer with a sense mass of a few kilograms and a baseline of X m 
can have a resolution of 10~ 4 E & er a 1 second sampling time (as 
reported in Wells [1984], p. 47). The primary error source for this 
sampling period is due to thermally induced vibration which is itself 
kept at a very low level since the instrument is operated at cryogenic 
temperatures. 



Figure 6. Canonical gravity grad ometer 
(from Wells [1984], p.43) 
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The Smithsonian Aatrophyaical Obaervitory (SAO)/ Piano Spaziale 
Nazionale (PSN) deaigni utilize gravity radiation antenna technology 
developed at PSN and cryogenic microwave cavity technology at SAO. 
The two approaches are being taken to determine which technology will 
be most suitable in a gravity gradiometer. Two single axis sensors are 
being constructed which will commence with the ultimate sensitivity 
goals of 10"* to 10"* E for the capacitive probe/radiation antenna 
gradiometer and 10 -4 to 10~* E for the cryogenic cavity gradiometer. 

The gradiometer being developed by Sperry Defense Systems could 
be described as a cryogenic version, of the Hughes design since they 
are related at least in principle. The Sperry design utilizes a 

magnetically levitated balance arm within a cylindrical bearing. The 
bearing itself is levitated and rotated at a constant rate such that the 
rotational rat# is one half the natural frequency of the balance arm. 
The gradient signal is generated through the detection of inductance 
created by the rotating system. The detection is made by using a 
superconducting quantum interference device (SQUID) placed at the axis 
of rotation. A SQUID is a parametric device whose output voltage varies 
in response to an input flux. The cylindrical bearing responds only to 
angular accelerations which are transmitted to the balance arm (which is 
sensitive to the gradients as well). Thus a feed back loop can be 

established to null the effects of the angular accelerations by monitoring 
the bearing's angular position. In the orbiting instrument concept, the 
cylindrical bearing is replaced with a rotating sphere and an ensemble 
of three such spheres would be grouped in a noncollinear configuration 
to constitute the tensor gradiometer. Precise orientation and altitude 
determination are a requirement inherent in this design. The noise level 
goal set by the investigators at Sperry is 3xl0“ 4 E. 

Completing the list of active gradiometer developers is the design 
put forth by the University of Maryland. This design also uses SQUID 
technology to detect changes in the magnetic flux which is coupled to 

the displacement of a proof mass. Paik summarized the principles 
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involved very neatly in one sentence (Paik (1981]): 

"an acceleration Ag drives the proof mass to a displacement Ax, 
which is converted to a magnetic flux signal At by means of the 
superconducting inductive transducer and finally the magnetic flux 
is detected by the SQUID producing a voltage output AV." 

The accelerometers can be placed in a configuration similar to that 
shown in Figure 6. When this configuration is used then the gradients 
are determined by differencing the outputs of the individual 
gradiometera following the concept underlying the definition of an 
Eotvos. In Paik [1981], ways of determining these differences are 
discussed and the concept of current differencing is introduced. 
Current differencing techniques provide a conceptually simple and 
efficient method to determine the in-line and cross gravity gradients. 
In the Univ. of Maryland design, the Niobium proof mass is weakly 
suspended which will limit gradiometry sensitivity to about 10~ 4 E over 
a 3 second integration time. Paik notes that the free-mass scheme 
would allow higher sensitivities and is realizable in space. Problems 
arise in attempting to test such a design in a "hostile terrestrial 
environment", i.e. the laboritory. Paik indicates further that a 
free-mass design in the zero-g evnironment of space could provide 
sensitivities up to 10“* E, 

In conclusion, it should be noted that cryogenic gradiometera 
promise higher resolutions than the conventional designs. However, the 
conventional designs for terrestrial and airborne applications are in a 
more advanced state of development. It should be noted too, though, 
that the latest push in gradiometric instrumentation is towards 
cryogenic designs and given enough time, with appropriate funding, 
cryogenic gradiometers could be available, in the next few years, for 
Earth gravity field determination as well as gravity missions to other 
terrestrial planets. This, then leads into a discussion of mission 
characteristics and proposed time tables. 


2,3 Gradiometer Mission Characteristics 


In this section, some aspects of mission and spacecraft design will 
be addressed. Many options exist at the moment and specific decisions 
are scheduled to be made over the next few years. Thus, only 
recommendations and assumptions will be made to delineate mission 
parameter limits to be used in the remainder of this study. 

Since the principle objective of the gravity gradiometer mission is 
to obtain a global data set of gradient tensor values, the satellite should 
thus be placed in an orbit to facilitate global coverage, This is most 
easily done by placing the satellite in a polar orbit with the inclination 
angle to be as near to 90* as possible. Furthermore, as was mentioned 
in the introduction, higher resolution of quantities of interest (e.g. 
gravity anomalies) is acheived with lower altitudes. The atmospheric 
drag experienced by the satellite at these low altitudes becomes the 
primary factor for premature termination of the mission. Two 
approaches have been proposed. The first, as mentioned earlier, is the 
use of the Tethered Satellite System in which the instrument package 
with the gradiometer is lowered from the command vehicle (in this case 
the space shuttle) and measurements are taken. Three main 
disadvantages arize with this concept which limit its usefulness strictly 
to the gradiometer test phase. These are due to the inclination of the 
orbit, and duration of measurmentB, and untested tether dynamic models. 
The space shuttle normally operates in a relatively low inclination orbit 
(30* -SO*) so that global coverage is impossible. The measurement 
density will likely be globally non-uniform due to the shortness of 
typical shuttle missions, which is undesirable from an analytical point of 
view. Finally, although investigations into the dynamic behavior of the 
tether systems are being made, the resulting models are largely 
untested in the real space environment. Until such tests are made, the 
gradiometer system's response to the tether's dynamic behavior remains 
uncertain. 


The second aproach to overcoming the effect of the atmosphere for 
low orbits borrows technology from the GRM mission. In order to 
maintain an orbit at (say) 200 km altitude, it is necessary to compensate 
for the dra»’ by employing thruster rockets. With a properly designed 
gradiometer, the accelerations from the thrusters will not affect the 
gradient outputs. It is this fact which causes gradiometers to be of 
interest in other dynamic applications (e.g. inertial navigation sytems, 
airborne and terrestrial gravity surveying). This thruster aided 
concept is formally called DISCOS (Disturbance Compensation System). 
DISCOS, as used in the GRM mission, is under rather tight tolerances 
and constraints. For the gradiometry mission, the constraints can be 
eased somewhat, since it is not required that the gradiometry package 
be treated as a drag free proof mass, sis is required in the satellite to 
satellite doppler tracking concept of GRM. It has been recommended 
that the free flying GGM mission utilize DISCOS technology or a variation 
of DISCOS. See Wells (1984), pages 59-60 for a more complete 
discussion. 

Data coverage is a function of satellite inclination and altitude as 
well as of time. For the best distribution of data it would be advisable 
to avoid altitudes which yield resonant orbitB with short period repeat 
cycles. The quality of distribution improves with longer mission 
durations. The optimal data distribution depends on the scientific goals 
of the mission. To guarantee that mission goals be met, it is advisable 
to overdetermine the system solutions (i.e. have too much data). 
However, mission planners must decide how much overdetermination is 
enough, given the economics of the mission. Also, in some cases, mission 
goals are relaxed. In this case it is easy to filter or smooth the raw 
data to meet these needs. 

In Figure 7, the effects of varying mission durations can be seen. 
The data was generated with the same program as was used in the 
construction of Figures la and lb. What it illustrates is that the 
recoverable X*xX* mean gravity anomaly accuracy improves as the 
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duration increases (i.e. more data is available). It ia quite easily seen 
that the curve levels off near 6 months and that a mission of twelve 
months only improves the mean gravity anomaly accuracy by 
approximately 6% in comparison with a six month mission. Thus it 
should be noted by system planners that extended missions may not 
realize significant improvements in determining static gravity anomalies. 

The demands for satellite ephemeris and altitude monitoring will 
provide the designers of the autonomous mission some technical 
challenges. For example, the uncertainty in the attitude sensing system 
must not exceed 3x10“ 7 radians in order to prevent degradation of 
measurement precision. In constrast, the space telescope pointing 
requirement is an order of magnitude smaller, at less than 3xl0~* 
radians over a twenty minute observation period. It is expected that 
the attitude control and monitor problems can be overcome as mentioned 
in Wells [1984]. The design put forward by Paik (1981] should be 
capable of monitoring attitude and position directly by integrating 
certain gradiometer and accelerometer outputs. These integrated 
quantities can further be placed in an absolute reference frame by 
using state of the art conventional attitude sensors ana tracking 
techniques. 

It is necessary that certain assumptions be made concerning the 
parameters associated with an automomoue satellite gravity gradiometry 
mission for the purpose of the primary investigation contained herein. 
These assumptions are outlined in Table 2. Using these parameters, 
then the raw data should be distributed every 4' in latitude and every 
3.6' in longitude at the equator. The east-west spacing of the data 
decreases towards the poles due to the convergence of the orbits at the 
poles. The sensitivity of the instrument (iO“*E) was selected as a 
conservative estimate. Up to this point, the sensitivity goal of 10~ 4 E 
for the mission has not actually been achieved by any of the 
organizations developing gradiometers. The investigation scientists 
remain optimistic however, that sensitivity goals will ultimately be met. 


Table 2 , 

Assumptions for Mission Parameters and Characteristics 

* Mission parameters 

* Orbit - low eccentricity (circular) 

* Inclination near 90 (polar) 

- altitudes of 200 km and ICO km 

* Duration - six months 

* Sample rate - i set of tensor measurements per second 

* Instrument assumptions 

* Cryogenic full tensor grttdlometer of the University of Maryland type 

* Sensitivity - 10*3 £ pg a ] ~ 1 ) 

* Attitude and ephemeris assumptions 

* Data assumed to be corrected for attitude and attitude rates and 
converted from electrical output signals to real physical values 

* Data assumed to be rotated Into local level cartesian coordinate system 

* Data assumed to be geographically tagged (latitude, longitude, altitude) 


The particulars concerning the attitude and ephemerie assumptions are 
beyond the scope of this work but these assumptions will need closer 
examination after specific design decisions have been made. The 
rotation of the gravity gradients, from the orientation at the time of 
measurment to the local level coordinate system can be made by time 
tagging the gradient output signal and the absolute attitude information. 
The local vertical can be determined by the gradiometer output and the 
required rotation angles can be determined from the attitude 
information. The geographical tagging can be made by several methods 
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including integration of gradiometer outputs, supplemented with global 
positioning system (GPS) simultaneous phase measurements (described by 
Melbourne t Tapley [1983]) and other conventional tracking methods 
including laser and unified *-bcnd radar tracking. 

2.4 Previous Reduction Schemes 

Several investigations have been made concerning techniques to 
reduce the satellite altitude gravity gradient measurements to some 
geophysically useful form. Two basic classes of reductions can be 

defined; global solutions and local solutions. The proposed global 
solutions attempt to determine the coefficients of the spherical harmonic 
expansion of the disturbing potential. The proposed local solutions 
attempt to determine directly, gravity quantities such as mean gravity 
anomalies and mean height anomalies (or mean geoid heights). A brief 
synopsis, highlighting certain aspects of these previously proposed 
solutions, follows. 

Many gradiometer studies were made in the 1960’s, however, most 
of these were proof of concept studies and dealt primarily with the 
improvement of inertial navigation system error sources. Geodetically 
relevant gradiometer studies began formal investigation in the early 
197f)’s. Since then, the discussions have followed two paths; airborne 
gradiometry and spaceborne (satellite) gradiometry. Much work has 
been done for the airborne application and recent investigations by 
Jekeli [1985] illustrate the immediate usefulness of the ur borne case. 
Although the airborne application and associated reductions will not be 
included here, it is important to recognize the close linkage between the 
two cases. The primary differences lie in the types of usable 
instrumentation and in the height of the gradient measurements. 

One of the earlier investigations of satellite gradiometer data 
reduction appeared in a thesis by Glaser [1972]. In this work, Glaser 
proposed an algorithm which computes improved estimates of the 


coefficient* in the spherical harmonic expansion of the earth's 
gravitational field. Hence, Glaser's solution is of a global nature. It is 
very much a preliminary work as it was not possible for Glaser to fully 
test the algorithm. However, the proposed technique was novel for its 
time. Essentially, Glaser proposed that the measured gradients (from a 
rotating instrument) be downward continued to the earth's mean radius 
and then integrated to directly yield improved coefficients. Probably 
the factor most likely to have hindered Glaser's attempt at programming 
the algorithm was the considerable computational effort required. Even 
with the prospect of computational difficulty, Glaser's proposal, though 
very bold, set the starting point for further investigations. 

That same year (1972), Reed outlined a reduction scheme which is 
less demanding with respect to the computations (Reed [1972]). Reed 
proposed a local solution where classical least squares techniques are 
used to determine 2*x2* and 5*x5* mean gravity anomalies referred to a 
14x14 degree and order geopotential model. The results presented are 
in the form of an absolute error analysis, where the 2*x2* (and 5*x5*) 
mean anomalies predicted from the simulated gradients at altitude were 
differenced from the gravity anomalies used to create the simulated 
gradients. Reed was very perceptive in making investigations for both 
strap-down and rotating types of gradiometers, although his conclusions 
concerning the future of strap-down systems was faulty. Much of the 
present report draws from the analytical work of Reed. 

The next definitive study on the geodetic usefulness of satellite 
gradiometry was made by Krynski and Schwarz [1977]. They were 
among the first to utilize the technique of collocation to compute 
geodetically relevant information from satellite gradiometry. In their 
study, they did work on determining the appropriate covariance model 
parameters, configuration of the gradiometer signals, and on combined 
solutions using additionally, terrestrial data. Their results, primarily in 
the form of a very generalized error analysis, yielded some fairly 
interesting conclusions. Among the more interesting of these is the 
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conclusion that only those gradients containing at least one radial 
derivative contribute significantly in recovering gravity information (in 
this case geoid heights). That is, only the T,„ T,„ and T„ gradients 
in the local level system contribute significantly. What appears to be 
missing from their study, however, is a discussion on the 
non-stationarity and anisotropic nature of the gradient covariances. 
Their conclusions remain sound in the case that their adopted 
covariance function behaves realistically. 

Recently, a new attempt at the global solution has been outlined by 
Rummel & Colombo [1985]. Due to the remarkable breakthroughs in the 
computational methods of harmonic analysis on the sphere by Colombo, a 
new interest in large global solutions has resulted. Their propoeed 
solution solves for the coefficients of the spherical harmonic expansion 
in an iterative process which includes orbit displacement and orientation 
uncertainties. A prerequiste for using the advanced analysis algorithm 
(decribed in elegant detail in Colombo ( 1981 ]), is that the measured data 
be in gridded form. Rummel & Colombo suggest creating cells of data 
averages which are then ’’dropped" onto a surface of revolution which 
better aproximates the actual orbits to create the gridded data set. By 
placing the data on an orbit approximating surface (rather than a 
sphere) helps to decrease the number of iterations required due to the 
reduction in the initial orbit displacement error. Rummel & Colombo ran 
computer simulations which showed, quite positively, the success of the 
method. There is much promise for this method because it has overcome 
many obstacles, notably the reduced computer time for the solution. 
The computer time required for the global solution is expected to remain 
substantial (i.e. greater than a few hours) thus, local solutions for 
direct gravity parameter recovery, such as the one outlined in this 
work, remain useful for those interested in regional investigations with 
a more direct approach (i.e. gravity gradients ■* gravity parameters 
rather than the three step method: gravity gradients ■* geopotential 

coefficients •+ gravity parameters). 








3. APPLICATION OP LEAST SQUARES COLLOCATION 


3.1 The Collocation Model 

The techinque of least squares collocation has been used 
successfully for a wide variety of problems in physical geodesy. It has 
been especially useful in certain satellite mission applications uuch as in 
satellite altimetry (Papp [1963]), satellite-to-satellite tracking (Hajela 
[1981]), and in estimating gravity potential differences between 
continents from satellite laser ranging data (Hajela [1983]). The results 
from these reports showed that significant solutions can be found with 
limited observational data. As is well known, for large observational 
data sets, the collocation technique suffers from the requirement that a 
large matrix be inverted or a large system of equations be solved. This 
drawback can be bypassed in two ways. First, the simpler approach is 
to determine what elements of the total data set can be ignored or 
which elements can be smoothed or averaged such that the errors 
resulting from the collocation aoluton remain small. In this approach, 
one may be forced to accept some loss of resolution in order to reduce 
the computation time. The other method is to switch to frequency 
domain collocation where the inversion con be handled with less 
computational burden. Details of the latter technique are found in Eren 
[1980]. The present study will follow the former, more simple approach. 

The fundamental equation of the least squares collocation technique 
is (Moritz, [1980], p. 102) 
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where 4 is ths vscior of observations, C 44 = (Ctt ♦ C nn ) is the sum of 
ths covariances of ths observations and ths observational noise, C B t is 
the covariance otatrix relating the observed quantities to the quantities 
being predicted, S is the vector of predicted quantities. 

The observational noise is usually associated with the observed signal 
uncertainty and in most cases, this uncertainty is treated as 
uncorrelated from one observation to the next, thus the C na matrix will 
be diagonal with all the elements on the diagonal having similar values. 

The collocation technique is especially appealing to geodesists since 
it allows for an estimate of the inherent error convariance matrix 
associated with the estimate vector S. The error covariance matrix (E aa ) 
is found through (Moritz [1980], p. 105) 


Bii = c ss “ C at C 4 ,Cl t 


(III. 2) 


where all quantites were defined above and, in addition, C aa is the 
covariance matrix of the estimated signal. The diagonal of Ef § provides 
the variance of the estimated quantities S. In the case that only one 
signal is being predicted (i.e. S has only one element), then Egft will 
ha*'e only one element and this value will be the variance associated 
with S. 

Another aspect that makes collocation a desirable technique is due 
to the fact that various data types can be combined to yield perhaps 
more improve estimates of the quantities desired. This aspect will not 
be a factor in this study because the intention of this study is to 
examine one specific data type (i.e. gradiomctry data). 


3.2 Bxplictt Matrix Structure-Observation Vector 


In the following discussion, the developments will be made in the 
most general way possible. This means that the "general case" 

corresponds to the case that all five independent gradient tensor 
components are uaed in the solution. If fewer components are used, as 
is the case in this investigation, it is not difficult to modify the 
following development to accomodate the changes. 

The vector of observations, I will then be composed, in the general 
case, of the five independent residual tensor component values 
determined at m particular location. To obtain these residual values, 
several preprocessing steps are required. 

First, it is assumed that the instrument outputs have been 
converted to physical units and that these values have been rotated 
into the local level coordinate system. The angles required for the 
rotation are determined from the satellite attitude control subsystem. It 
is further assumed that the attitude information has a propagated noise 
level below that of which is detectable by the gradiometer, thus the 
rotation angles are assumed errorless. Once these preprocessing steps 
have been completed, then what results are the measured gravity 
gradients, uhown as a "vector segment" below, 

VjS 

V 

v” (III. 3) 

Vis 

V ” J {? i. Ai, ri } 


where the numerical subscripts denote the differentiation axes 
corresponding to Figure 4 and the quantities in braces are simply 
reminders that the measurements are position tagged with Aj, rj 
being latitude, longitude, and radial distance from the geocenter for the 
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ith observation respectively. 

Grouping of the vector segments can be made in many ways, but for 
the purpose of gravity parameter estimation, the logical alternative ia to 
group by geographic location. If the data is randomly distributed, it 
car be grouped by choosing the spherical distance f 0 and including aa 
segments in the observation vector, those gradient values at 4j, Xj such 
that 


cos ‘[sintp sintt + cosip cos ( X i - X p ) ] * y 0 


(HI. 4) 


where ip and Xp are the latitude and longitude of the prediction point. 
However, using randomly distributed gravity gradient data causes 
problems in computing the required covariances. One way to overcome 
these problems is to grid the data in some fashion. 

There are several data gridding schemes available and depending on 
the type of solution (i.e. global or local) some schemes may be 
excessively complex or simplistic for application here. The first step in 
gridding the data is to predict point gradient values at the grid 
intersections of a regional block (cf. Figure 12). A least squares 
collocation predictor could be used for this task where nearby 
gradiometer measurements are used as input to form the predicted 
values. Since the gravity gradients, as will be shown later, are a short 
wavelength phenomenon, the prediction should be based upon only those 
measurements that are sufficiently close enough to the predicted grid 
value to keep the loss of gravity information to a minimum. The grid 
values should be reduced to a specific surface in order to ease the 
covariance computational demand® in the main collocation solution for 
gravimetric quantities. Moving the height dependent aspect of the 
gradient signal to the preprocessing step not only simplifies the 
computation of the covariances in the main collocation algorithm, but also 


provides the possibility of forming the observation vector as will be 
described shortly. The reduction to the specific surface can be 
accomodated in the prediction step. The reducing surface assumed in 
this study is a sphere with a radius of R N from the geocenter, which is 
equal to the average goocentric radii of the gradiometer measurements 
within the regional block. Thia surface should be sufficient in our case 
since the orbit is considered circular and we are considering 
observations in a relatively small region or block (with, in general, Vo 4 
5*). Another auch surface is a surface of revolution about the earth's 
axis which better fits the actual orbits (Rummel k Colombo (1985]). This 
surface is not used in this work because the radial uncertainties which 
are being solved for in Rummel I Colombo’s global solution can be 
ignored in this type of solution without detracting from the overall 
results. It is also not used because of the complications resulting from 
the height dependent covariance computations. Thus, once these 
averages and reductions are made, then what results is a geographical 
grid of points, to each of which are attached the five predicted and 
reduced independent gradients shown in Figure 8. 

The complete gradient vector is then structured from point 1 to 
point MN. In transposed form the gradient vector is 


< T = tVi s , Vja, Via, V}*, V$£, V?S, VfJ, V?S, V?5f] , (III. 5) 


where the superscripts indicate which point, the gradient value pertains 
to. Thus the observation vector, for the case that all five independent 
tensor components are entering into the solution, will have 5(MN) 
elements. For the formulation used by Krynski k Schwarz [1977], only 
the three components containing at least one radial derivative enter into 
the solution. Thus their gradient vector had 3(MN) elements of the form 


j V * 
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* T = [Vi,* 


V ! 


ts * 


vi vSff, v?S, VS}] 


(III. 6) 


For the usual case that the region considered is equi-angular, then 
At will equal AX and M=N (in Figure 8). Equi-area regions should be 
used in the polar regions due to problems associated with the 
convergence of the meridians. 
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Figure 8. Numbering scheme for gridded sub-blocks 


To complete the observation vector set-up, the differences between 
the predicted grid gradient values and the computed reference 
gradients are taken to form the residual gradients (or anomalous 
gradients). This is done to meet the requirement that the observed 
data be centered. To center the data, reference point gravity gradients 
are computed for all components being considered in the solution at 
each point (4|, Xj) on the grid at radius R&f. The reference gradients 
are computed from a aet of reference geopotential coefficients (C^ m > 
Snm) complete up to degree and order N. The reference coefficients, 
together with (if, X{, R N ) are inserted into the properly differentiated 
version of equation (II.4). The differentiations were made by Reed 

[1972] and are summarised in Appendix A. Since the required reference 
gradients are on a grid, then advantage can be taken of the 
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mathematical symmetries which result, and a computationally efficient 
algorithm can be utilised, auch aa described by Robbinw [1984]. The 
observation vector, ♦ can then be establiahed by (cf. eq. (II. 5)) 


|i = 


Vis 

Vis 

Vis 

Vl, 

Vi, 


Ui, 

Uis 

Ui, 

Ul, 

Ui, 



Ti, 


Tit 

II 

Tit 


Tl, 


Tis 




where I* is the observation vector segment. 


(III. 7) 


3.3 Explicit Matrix Structure - Covariance Matrices 

The overall structure of the covariance matrices is dependent upon 
how the observation vector and the resulting signal matrix are 
structured. The basic structure of the observation covariance matrix 
will be 5MNx5MN, where MN are the total number of points contained in 
the NxM grid of locations. The matrix illustrated in Figure 9 will need 
to be inverted according to equations (III. 1 ) and (III.2) . Since the 
matrix is symmetric, some inversion computation time can be saved by 
storing the matrix in symmetric storage (vector) mode. However, for the 
purposes of testing and keeping computation times relatively short (lesB 
than 10 minutes on the IBM 3081 at the Ohio State University), the 
author selected 1000x1000 as the maximum dimension of the 
autocovariance matrix. If all five of the components are used (as shown 
in Figure 9), then, at a maximum, a grid of 200 points can be 
considered. This corresponds to a 14x14 grid of residual gravity 
gradient observations. More points on the grid can be utilized when 
fewer than five components enter into the solution. 

The covariance matrix relating the observed quantities to the 
desired signal, C B t» has its structure defined through the number of 
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observations and the number of predicted quantities. It will have the 
dimension L x 5MN where L is the number of predicted quantities and 
5MN is the total number of observed gradients in the grid (general 
case). The predicted quantity can be any gravimetric quantity (e.g. 
gravity anomalies or height anomalies, as used in this study). 

In general, C B t will be a row vector because normally only one 
quantity is being predicted. That is, only one gravimetric quantity, 
located in the center of the grid (and on the earth's approximating 
sphere), is computed for a specific grid. For adjacent predicted 
gravimetric quantities, the grid is shifted and observations symmetric to 
the computation point enter into the computation. Thus, a moving 
window technique is used in which the covariances are longitudinally 
invariant. Computation of the covariances themselves, follows. 


3.4 Covariance Computations 

The covariances associated with the gravity gradients as well as the 
gravimetric quantities can be expressed as linear functionals of the 
disturbing potential. This well known property contributes to the 
success of collocation. The general form for the covariances of the 
disturbing potential, T is 

I N 

IZ H£T (2n+1) £ fl p n( c ° 8 fpp') 

n - 2 

m 

+ L_ H£') n (2n+l) <rflPn(cos -*pp' ) } (III. 8) 
n=N+l rr ' 


where r, r ' are the radial distances to point P and Q respectively, yp p ' 
is the spherical angle separating P and Q, a is the mean equatorial 
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•re the Legendre polynomials of degree n. 


£* = r (* a Cn» + **S|ni)/(2n ♦ 1) (III. 9) 

»=o 


with the absolute error of the Modeled coefficients given by 


cCng = CnB(tru«) ~ C UM(mod«l) (III. 10a) 

kSqb = S nB ( iru# ) - S nB ( Moda |) (III. 10b) 

Since cC nD) and cS nm are unknown, then aa an approximation, the 
coefficient standard deviations of a spherical harmonic geopotential 
model solution, eC nm and eS nm are substitued for eC nm and eS ntn . The 

corresponding term, efl in the right hand side of (III.8) is the degree 

variance of the potential coefficients; 


= ZI (Cfl. + Sfl»)/(2n + 1) 

M=0 


(III. 11) 


This formulation is originally due to Colombo [1980] (section 3.2) but the 
notation of Hajela [1983] (also section 3.2) is used here. The potential 
degree variances are related to the anomaly degree variances by (Hajela 
[1983], p.17): 
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c n = rS(n - l) a (2n + l)<r», (III. 12) 

where y 0 is a mean value of gravity. 

The anomaly degree variances, c n , are usually given by empirically 
derived models* Many collocation studies have used the Tscherning & 
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n (n-2) (n+B) 


(III. 13) 


where a and B are model parametera. The main reason for its extensive 
use is due to the availability of advanced software baaed upon this 
model (see Tscherning k Rapp [1974], Tscherning [1976), and Sunkel 
[1979]). This model exhibits logarithmic behavior and as such, it is not 
the most suitable model for certain types of covariance modeling. It has 
been especially criticised for its unreasonably high gradient variance 
especially in association with the parameters determined by Tscherning 
k Rapp (1974). These parameters were determined by least-squares 
fitting the model, equation (III. 13), to observed anomaly degree 
variances from degree 3 to 20, 1* and 6* mean (block) anomaly 
variances, and a point variance of 1795 mgal 2 . More recently, attempts 
to update the parameters were made by Rapp [1979], and Jekeli [1978] 
attempted using various other parameters to determine their effects on 
the variances of several gravimetric quantities (including the vertical 
gravity gradient). The specific parameters and the effects upon 
gradiometry solutions will be discussed in the next chapter. 

Other models for the anomaly degree variance have been proposed 
by several investigators, we mention only the two component-model 
attributed to Moritz. It's basic form ic 


(III. 14) 


where the <x,, a a , S,, S 3 , A, and B terms are the model parameters. 
Parameter investigations were made by Jekeli [1978], Hein and 
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Jochemczyk [ 1979], and Rapp (1979). Although it haa been ahown by 
theae investigators that more natural anomaly and gradient variances 
can be obtained from the two component model, this study will restrict 
itself, with regard to the analysis, to the Tscherning and Rapp model 
for the aff or mentioned reason of readily available advanced software. 
The results contained herein might be improved by using the 
two-component model. 

The particular covariances are derived through the law of 
covariance propagation. The gravimetric quantities, i.e. gravity 
anomalies and geoid heights are related to the disturbing potential by 
(Moritz (1980), p.108) 


4g = - 
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(111. 15) 

(111. 16) 


with ro being a mean value of gravity. Both equations are given in 
spherical approximation. The functionals relating the gradients in a 
local level coordinate system to the disturbing potential are slightly 
more complicated. The functionals for the tensor diagonals are (Reed 
(1972), p.32); 


T *> = Fw? Tax - *3* T» + f T r (III. 17a) 
T aa = pr T** + ~ T r (III. 17b) 
Taa = T rr (III. 17c) 


The sum of the above expressions yields Laplace's equation in spherical 


coordinates. The off diagonal tensor functionals are: 


T., * T„ * r , c ^ T X , * T» 

(III. 18a) 

T,s " T »* ’ rcost Txr " r’cost T * 

(III. 18b) 

T„ « T„ = £ T» r - ^ T* 

(111.18c) 


Now applying the law of covariance propagation with Lj and Lj denoting 
the functionals, the covariance of particular interest is (Moritz [I960], 
P»87) 


Hij(P, 0) = LfLf'K(P.Q) 


(III. 19) 


where K(P, Q) is given by (III.8). Substituting the functionals in 
(III.15) through (III.18) into (III.19) will yield all of the covariance types 
required in the "general case solution”. These expressions have been 
derived by the author and are given in Appendix B. 

It should be emphasized however, that in some cases, the resulting 
covariance expressions are no longer isotropic. That is to say, they are 
no longer strictly functions of the radial distances, r and r', and the 
spherical distance, Tpq* The gradient covariances can be treated in a 
manner similar to the treatment of vertical deflection covariances 
described by Tscherning k Rapp [1974]. Recently, Krarup k Tscherning 
[1984] have recast the T, 3 , T aa and Ti a gradient component covariances 
into isotropic form for use in collocation solutions utilizing torsion 
balance observations. Tscherning has further modified his closed form 
covariance algorithm (COVAX) to impliment the isotropic form (Tscherning 
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(1963)). 


3.6 Reduced Forms of the Collocation Solution 

Since Krjrnski I Schwarz [1977] concluded that the gradients 
contributing significantly to the local collocation solution are thoae 
containing at least one radial derivative, then the "general case" of 
utilizing all five of the observable gradients can be reduced to utilizing 
the three radially differentiated gradients (T,,, T, ,, T aJ ) with negligible 
loss of accuracy. If the same number of grid points is retained as was 
used in the consideration of the "general case" solution,) then the 
computational burden and necessary computer time will be reduced 
especially for the in verson of C t$. In this case, the observation vector 
will be 



vi. 

- ui. 


Ti, 


«i = 

vl. 

- u}. 

a 

Tj, 

(iei*MN) 


Vi, 

- Ui, 


Ti, 



(III. 20) 


where indicates the observation vector segment with i being a 
number associated with a set of observation on a grid and where MN is 
the total number of points in the grid ( ct . Figure 8). The observation 
covariance matrix will be similar to that shown in Figure 9 but without 
the T aa and T ta related covariances. The dimension of Cn is, in this 
case, 3MNx3MN and if it is again assumed that the maximum dimension of 
practical inversion for the purpose of testing is 1000x1000, then the 
grid can be enlarged to the size of 18x18. The enlargement can be 
made by densifying the grid (of specific geographic coverage) or by 
retaining the data density and extending the geographical coverage. 
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The formulation and covariance computations are simplified further 
by only considering the second radial derivative component (the vertical 
gravity gradient, T>>). The covariance function for the vertical gravity 
gradient is isotropic from the start sinco the radial derivative 
functional, when applied to the anosialous potential covariance through 
the covariance propogation law, does not disturb the isotropic nature of 
the anomalous potential covariance function, K(P, Q) (cf. equation 
(III.19)). Only one observational type needs to be considered therefore 
only one autocovariance table needs to be generated. The observation 
covariance matrix will have dimension MNxMN and with the assumed 
maximum dimension of 1000x1000, then the grid can be further enlarged 
to 31x31. Again, either the density of observations can be increased (to 
a limit not to exceed the actual measurement distribution) or the 
geographical coverage can be made more extensive than for the case 
mentioned in the previous paragraph. 




4. RESULTS OF THE ERROR ANALYSIS 


4,1 Description of the Investigations 

The investigations reported herein are primarily in the form of an 
error analysis* This is due to the fact that at this time, a satellite 
gravity gradiometer mission has not yet been tested, so no real data yet 
exists* The basis for the error analysis comes from the expression of 
the error covariances (equation (III. 2)) which is repeated here for 
reference (ci. section 3.1 for term definitions), 


= C»b "CitCiWt 


(IV.l) 


This way the errors of the predicted signals are found by the model 
implied covariances, which are used to form the matrices C BB , C B t, and 
Ctt, and the instrument noise level which enters through C nn (recall 
that C|| s Ctt * ®nn). 

The original idea was to impliment the "torsion balance version" of 
Tscherning’s COVAX program to compute the covariances for those 
gradients containing at least one radial derivative. However, the results 
of this attempt proved unsatisfactory either due to the instability of the 
inversion of C<| or due to undetected software errors. Thus to 

simplifiy the investigation, it was decided that only gradient components 
with two radial derivatives would be used as the assumed observable 
(referred to as the radial-radial component, T ss , in the remaining 
discussion). 
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The remainder of this chapter will include a diecueeion on the 
generation and behavior of the covariance*, and presentations of the 
error analysis with variations in the parameters. 


4.2 Covariance Generation and Behavior 

Several methods have been devised to generate covariance tables. 
Two of the more popular methods are by truncated series and by closed 
form expressions. The truncated series method has been used 

occasionally with the two component degree variance model of Moritz 

(equation (III. 14)). The closed form expression method has been 

programmed by Tscherning and Rapp [1974] and further modified by 
Tacherning [1976] and [1983]. The degree variance model used in the 
closed form expression has traditionally been that Tscherning & Rapp 
(equation (111.13)), however recently the Moritz model has been 

implimented in a closed form algorithm by Hein [1981]. 

The closed form algorithm of Tscherning [1983] has been used 
throughout the present work. This algorithm (referred to as COVAX in 
thin report) was selected since it was readily available to the author 
and since it is easily adaptable for the use of higher order reference 
potential models. 

The selection of parameters for use in the Tscherning & Rapp 
degree variance model poses an interesting problem. The parameters 
computed by Tscherning & Rapp [1974] implied a horizontal gradient 
variance of 3500E 3 which is considered unrealistically high. The gravity 
anomaly variance of 1795 mgal a implied by their parameters has been 
accepted for several years but recently the value has come under close 
scrutiny. Evidence from studies of gravimetric quantities and empirical 
covariances in Canada have suggested that the 1795 mgal a value may be 
too high. For the Canadian region, Schwarz & Lachapelle [1980] 
estimated the gravity anomaly variance to be 837 mgal 3 . One must 
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remember though, that the Canadian evidence ia of a regional nature 
and thus it is limited to an incomplete set of gravity anomaly producing 
geophysical structures (i.e. the samples taken in Schwarz k Lachapelle’s 
investigation do not include a global set of geophysical structures such 
as deep ocean basins* trenches, island arcs etc.)* Thus their result for 
the gravity anomaly variance may be somewhat low, perhaps not by a 
very large amount. The horizontal gradient variance was estimated to 
be approximately 20GE* in Canada by Schwarz & Lachapelle which is in 
fair agreement with slightly higher gradient variances determined by 
Hein k Jochemzcyk [1978] in Germany. Therefore, the best degree 
variance model parameters to use are those which imply a global anomaly 
variance of approximately 1100 mgal 3 (Rapp, private communication) and 
a horizontal gradient variance of approximately 300 E a (as a compromise 
value between the results of Schwarz k Lachapelle and Hein k 
Jochemzcyk). One set of parameters for the Tscherning k Rapp degree 
variance model meets these variance requirements, which were 
determined by Jekeli [1978]. These parameters, as well as other used in 
previous investigations, are listed in Table 3 (cf. equations (III. 13) and 
(III. 14) for the expressions of the degree variance models). To illustrate 
the power spectrum of the three parameter sets of Table 3, the products 
of the degree variances with their respective S-terms are plotted by 
degree in Figures 10o and 10b. This product will be called the scaled 
degree variances and are found by 


c n = c n sn+a 

In = g n S" +3 = c „S n+3 


(IV. 2) 
(IV. 3) 


with c n given by the Tscherning k Rapp model (equation (III.13)) and 
with g n denoting the vertical gradient degree variance. For the Moritz 
two component model, the scaled degree variances are 
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c n = c n ( c n from equation (III. 14)) 



«t §TT + 


n+A 


*1 


(n~0 

(n-*)(n+*) 



The variances are found through the summation 


(IV. 4) 
(IV. 5) 


Co = 


m 


(IV. 6) 


Cqv 


I In * 2 G 0 m 
n=a 


(IV. 7) 


with the last approximate equality shown by Jekeli [1978]. 

Figure 10a illustrates the scaled gravity anomaly degree variances 
for the three parameter sets of Table 3. The numbers associated to the 
curves correspond with the enumeration in Table 3. It can be seen that 
the Tscherning k Rapp model with Jekeli’s parameters exhibits lower 
power than the other two models for degrees below 800 this is the 
reason for the lower value of the gravity anomaly variance, especially in 
view of (IV.6). Also notable, is the similarity between the Jekeli 
parameter model and the Moritz with Hein parameter model beyond 
degree 800. Tscherning k Rapp's parameters exhibit greater power for 
the high degree spectrum. All curves can be seen to exhibit high 
power at degrees below 700. 

Figure 10b illustrates the scaled vertical gradient degree variances 
for the same parameter sets. The effects due to the two components of 


Table 


Uaed in thia work: Tscherning It Rapp aodel with Jekeli paraaeters 

a = 343.3408 agal 3 
B = 24 

S = 0.9988961 
Iaplied variances: 

(Co: gravity anaauily variance, Q 0H : horizontal gradient 

variance. ) 

Co = 1089.5 agal 3 
G oh * 338.9 B 3 

Tscherning & Rapp aodel with Tscherning It Rapp paraaeters: 
a = 425.28 agal 3 
B = 24 
S = 0.999617 
Iaplied variances: 

Co = 1795 agal 3 
Gqh = 3500 E 3 

Moritz two component aodel with Hein [1981] paraaeters: 
a, = 7.516922 agal 3 a, = 82.04054 agal 3 

S, = 0.994425 S, = 0.9996642 

A = -2 B -7 

Iaplied variances 
Co = 1800 agal 3 
G oh = 1000 E 3 







Figure 10a. Scaled gravity anomaly degree Figure 10b. Scaled vertical gravity gradient 

variances implied by the three parameter degree variances implied by the three para- 

sets of Table 3. 1: Jekeli parameters, meter sets of Table 3. The numbering schem< 

2: Tscherning & Rapp parameters, 3: Hein is the same as Figure 10a. 
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the Moritz model are strikingly evident by the two maximume reached at 
n=390 and n=3000. The rather sharp maximum at n=390 ia due to the 
logarithmic portion of the two component model. The two parameter eeta 
baaed upon the Tacherning & Rapp model display similar apectral 
behavior. It appears that the spectrum beyond (roughly) degree 500 
contributes to the high horizontal gradient variance implied by 
Tacherning k Rapp's parameters. Again, all curves exhibit the well 
known fact that gravity gradients are primarily a high degree (short 
wavelength) phenomenon. It should be remembered that the real world 
spectral behavior of the gravity gradient variance is still largely 
unknown. Perhaps with the influx of terrestrial and airborne 

gradiometry data, the parameters of the degree variance models can be 
refined further to incorporate this new information. 

The effect upon the degree variances from using a high degree 
and order reference potential model to remove reference gravity 
gradients from the measured gradients is illustrated in Figures 11a and 
lib. The DEC 81 reference potential model of Rapp [1981] has been 
selected for this purpose. Recall from Section 3.4 that the reference 
potential coefficient variances (e a S nm i e a C nm ) are substituted for the 
indeterminable absolute error of the reference potential coefficients. 
Thus equation (III.9) becomes 




n 

-JZ 


(e*C t 


a=o 


+ e a 8 m )/(Mn + 0 


(IV. 8) 


The error potential degree variances are related to the error anomaly 
degree variances by 


*-(9 ) a (n-l)(2n+l)ea 


(IV. 9) 


The scaled error anomaly degree variances are found by substituting £<. 
from (IV.9) for the C n in (IV.2) 


£ c = £ c S n+a (IV. 10) 

further, the scaled error vertical gradient degree variances are found 
by substituting £ c for C n in (IV.3) 



(IV.ll) 


Figures 11a and lib illustrates these scaled error degree variances 
implied by the DBC 81 coefficient degree variances. The DEC 81 degree 
variances are a modified set of degree variances not reported in Rapp 
[1981]. The degree variances used in this study were recomputed by 
(IV.8) up to and including degree 36. Beyond degree 36 the error 
potential degree variances were computed by the expression 


£fl = e J (C, S) n n>36 


(IV. 12) 


where e(C, S) n are approximate degree accuracy estimates given by 


elC Sin = c (^g)8 + 

e(C * S ' n 2ro(n-l)v/n + 


sampling nrrors 


(IV. 13) 


where e(Ag) is a global estimate of the accuracies of the l’xl* mean 
gravity anomalies used in determining the coefficients; for the modified 
degree variance computation, a value of 10 mgal was used; 9 is the 
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are gravity anomaly degree variances for the are vertical gravity gradient 

three parameter sets of Table 3. Numbering corresponding to Figure 10b. 

scheme corresponds to Figure 10a. 
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block siEe ( 1 *xl * * i.e. 0=1 # ), and y 0 ia a mean value of gravity. 

The ecaled degree variancea belonging to the three parameter sets 
of Table 3 are shown in Figures 11a and lib for n>180. It should be 
noted that the error degree variances abruptly alter the spectrum of 
the modeled degree variances. This is especially true for the anomaly 
degree variances of Figure 11a. 

As geopotential models improve, their coefficient accuracies will 
improve further supressing the error degree variance influents in the 
covariance computations. As a minimum case, one can consider the 
reference geopotential model to be errorless (i.e. the coefficients are 
equal to the true coefficients in equation (III. 10)). In this case (also 
considered in the next section), the degree variance spectrum has no 
power up to degree 180 and power beyond that degree as given by the 
empirical degree variance model. Considering the reference potential 
coefficients to be errorless is useful in determining the sensitivity of 
the resulting accuracy of the gravimetric quantities to the coefficient 
accuracies implied by the geopotential model actually used. A discussion 
of the results of this consideration is made in the next section. 

To summarize, the covariance model of Tscherning A Rapp will be 
used for the remaining investigations due to its ease of operation and 
availability. The parameters of Jekeli [1978] were chosen to be the most 
compatible with the latest estimates of the gravity anomaly and 
horizontal gradient variances. 


4.3 Results of the Error Analysis 

The error analysis primarily consists of determining the expected 
error of the resulting gravimetric quantity computed through the 
collocation technique for an assumed geographic area at altitude. A 
computer program has been written which implements equation (IV.l) by 
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living covariances computed from the "torsion balance version" of COVAX 
by Tscherning. The program is called GIFRAD, for which a listing may 
be found in Appendix C. 

There are three primary input variables for program GIPRAD, these 
being; the latitude, ♦; the overall grid size, D; and grid spacing, DELD. 
The latitude is necessary for the computation of the spherical distances, 
equation (B.3) in the appendix, since the grid is to be based upon the 
lines of latitude and longitude. The overall grid size and the grid 
spacing are illustrated in Figure 12. These parameters can be varied 
which then offers the possibility for investigating the behavior of the 
error estimate with regard to theae parameters. The gridded gradient 
values lie at the intersections of the grid lines, as shown in Figure 12, 
at a mean satellite altitude computed as described earlier in section 3.2. 
The gridded values of Figure 8 are treated as point values in Figure 12. 
The error analysis utilized an equi-angular grid, thus referring to 
Figure 8, AA=A* and M=N, therefore in Figure 12, N=M=\D/DELD)+1. The 
value of M and N are constrained to be odd thus causing the central 
point of the grid to be positioned directly above the computation point, 
P. This causes a symmetric distribution of data which can be utilized to 
decrease CPU time during program runs. 



Figure 12. Grid size and spacing definition. D is the grid width and 
DELD is the grid interval. 



The computed gravity an< j height anomaly accuracies are in terms of 
30' *30' means which correspond to a spatial resolution of approximately 
56 km. Since mean values were chosen, the covariances must reflect 
this choice. Several methods have been devised to compute mean 
gravity quantity covariances. A good theoretical discussion can be 
found in Jekeli [1978], whereby the Pellinen smoothing operator is 
applied. However, a considerably simpler approach has been advanced 
in Tscherning and Rapp [1974] whereby the Pellinen operator is 
approximated by a height dependent quantity^ Jekeli [1978] writes the 
mean gravity anomaly as 


Cov(Al Pt 4l Q ) =XZ 0A c n sn * ap n( co «y P o) 
n=a 


(IV. 14) 


Tscherning and Rapp [1974] write, in the same notation, the approximate 
formula as 


Cov(Ag p ,Alg) = YZ c n l (R^ f h)» ] n+a s n+ap n( coeVVc ) 


n=a 


(IV. 15) 


Thus, the Pellinen operator, Pft is approximated by 


« ‘ [(R?M 


n+a 


(IV. 16) 


Since the Pellinen operator is a function of the spherical distance f 0 , 
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fin 


1 


1-COSfo 


i: 


p n (t)dt 

coifo 


(IV. 17) 


where t=coaf 0 » then it refers to the enoothing within a cap of the size 
However, block aizea of 30'x30' are of interest, thua the 
corresponding cap size can be determined from 


* ■ PSA* 


O'*’* 


(IV. 18) 


where 8 is the size of the block. Therefore when 8 - 0.5*, then from 
(IV.18), io = 16' 56". 

The value of the height, h which best approximates the 0 n -function 
for i > o = 16' 56" is found by comparing the gravity anomaly variaiice 
computed from (IV.14) for the particular block size 0 to the variance 
computed from (IV. 5). Since, in the case of variance computations, P=Q, 
then ths* Legendre polynomial term becomes unity. Thus, the 
comparisons are based upon the following equations: 


N 


Var.Ui) =2Z fifcn* 1 ** 


(IV. 19) 


n=a 


»«■•<«> 'II 


(IV. 20) 


where N is the maximum degree of the summation. The comparison was 
made by programming equations (IV.19) and (IV.20) with N=2000 by a 
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program called COVRQUIV. For tha praaant investigation, the degree 
variances uaed were 


c n = Cc (from eq. (IV. 9)) (2*0*180) 

c n = Tacherning fc Rapp aodel with Jekeli parameters (n>180) 


The results of the comparison yielded a value for h of 6389 meters. 

The covariances required in the radial/radial collocation solution 
were generated by COVAXX where the error anomaly degree variances 
implied by the DEC 81 geopotential coefficients, the parameters of Jekeli, 
and the height h associated with the mean gravimetric quantites were 
provided as input. The covariances, in table form, were placed in 
program GIFRAD where too, the input parameters D and DELD (cf. Figure 
12) were specified. Another set of runs were made to test the effect of 
the geopotential model errors upon the resulting gravity anomaly and 
height anomaly accuracies. This was performed by means of omitting 
the coefficient errors in the covariance computations. This is equivalent 
to assuming a perfect geopotential field model. The variances of the 
predicted quantities (C BB in equation (III.2)) area as follows: DEC 81 

model errors included; 302.12 mgal a and 1.161 m a gravity anomaly and 
height anomaly variances respectively; perfect goe potential model to 
degree 180: 224.74 mgal a and 0.114 m a gravity anomaly and height 

anomaly variances respectively. 

The final results of both sets of runB are illustrated in Figures 
13a-13g. Using the above variances for the predicted quantities, it can 
be noted that the inclusion of gradiometer data causes a substantial 
improvement in the mean gravity anomaly and height anomaly accuracies. 
Four curves are associated with each figure with the curves having 
open box type symbols referring to gravity anomaly accuracies and 
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curves with croaa type symbols referring to height anomaly accuracies. 
Note that the respective accuracy scales are plotted on each side of the 
figures. 

Several observations can be made from the results. First, the 
accuracies improve with smaller grid intervals. This is not unusual 
since the number of geographical points increases as the interval, DELD 
decreases while the grid width D remains constant, that ia, more data in 
entering into the solution. The accuracies would reach their minimum, 
for a specified grid width, when the grid interval becomes infinitely 
small, or in other words, if the gradiometry data were continuous over 
the region. However, this minimum cannot be achieved for two reasons: 
one, the dimension of 0|« would become infinite, and two, continuous 
data is not realisable by current and planned gradiometers. 
Furthermore, the minimum grid interval is constrained by the mission 
parameters, most notable, the data acquisition rate, which determines the 
raw data geographical spacing, and by the assumptions concerning the 
gridding technique. 

Noteworthy too, is the rather large jump in th>v accuracies for grid 
intervals larger than one degree when the coefficient uncertainties are 
included. This is a natural outcome of the spatial resolution of the 
reference model which is approximately one degree. Therefore, the poor 
accuracies resulting from the geographically sparse data distribution is 
due to the influence of the coefficient uncertainites. Note that the 

"errorless" accuracies retain a smooth character for grid intervals 
beyond one degree of arc. 

The figures also indicate a general increase of accuracy as the grid 
width is increased for grid intervals less than one degree of arc. For 
instance, if one considers a specific grid interval (e.g. 30'), then as the 
grid width is increased two things are happening. First, more data is 
entering into the algorithm and second, the errors due to the omission 
of the region exterior to the grid japs being reduced. Both are thereby 
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Figure 13f. 30'x30' mean value (&g, f) accuracy estimatea for 3 *.5 grid 
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Figure 14# Relationship between gravity gradients and gravimetric 
quantities. Arrows point in direction of smoothing (damping 
of high frequencies) with spectral operators also shown. In 
the diagram, R{ is a mean earth radius, y 0 is a mean value 
of gravity at the earth's surface, and r a is the radius of 
the satellite at altitude. From Rummel (1975). 













effecting the accuracies, causing their improvement. It should be noted 
that the rate of accuracy improvement per increase in grid width can be 
seen in the figures to be decreasing as the grid widths are enlarged 
towards 2*. 

Additionally, it should be noted that, in general, the mean height 
anomaly accuracy does not improve as dramatically as the mean gravity 
anomaly accuracy considered as a function of decreasing grid interval 
spacings (i.e. introducing more data into the algorithm). This is easily 
attributable ot the relationships between the mean gravimetric data dn 
the gridded gradiometric data (Figure 14). The tranformation from 
gravity gradients to gravity anomalies is functionally equivalent to an 
integration whereby smoothing takes place according to the spectral 
operators given in the figure. The situation is similar for the 
transformation from gravity anomalies to height anonalies. The 
transformation from gravity gradients is, in effect, undergoing two 
smoothings with a combined spectral operator of 


B g 

(n+1) (n+2) 


(IV. 21) 


which behaves like l/n a . Since in Figure 13a-13g, shorter wavelength 
data is entering into the solution when the data density is increased 
(by reducing the grid interval), then, from the spectral operators, it 
follows that the height anomaly accuracies behave smoother than the 
gravity anomaly accuracies since the operator (IV.21) suppresses the 
higher frequencies to a greater extent than the gradient-to-gravity 
anomaly operator 

R E (n-l) 

(n.l)(n.2) < IV - 22 > 
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The accuracy eaiimatea were computed for several other altitudes as 
* given in Table 4» The overall improvement of the gravity anomaly 

accuracies is due to the strengthening of the signal in the medium 
wavelength gravity inforsiation. The accuracy will reach a minimum at a 
certain altitude which is a function of instrument sensitivity, data 
distribution and data coverage. The accuracies of the height anomaly 
changed only a negligible amount since the improvement of the gradient 
signal at lower altitudes occurs at frequencies higher than those which 
contribute the majority of the total height anomaly signal. 


Table 4. 30 '*30' mean gravity anomaly accuracies computed by the local 

collocation algorithm GIFRAD for a grid width of 2 %* and a data spacing 
of IS'. Coefficient uncertainties included. Instrument sensitivity: 
10"* B at 40* latitude. 


Altitude 

Af accuracy 

(km) 

(agal) 

200 

9.36 

180 

8.60 

160 

7.61 

140 

6.96 


Finally, the sensitivity of the accuracies to the coefficient errors 
can be clearly seen in Figures 13a-13g. Since the height anomalies are 
most affected by long wavelength information, they show marked 
improvement when the long wavelength error in the coefficients are 
removed. On the other hand, the gravity anomalies show only moderate 
improvement when the potential model uncertainties are removed. This 
behavior clearly indicates that a highly accurate geopotential model 
when used in the local collocation algorithm can significantly improve 
the resulting gravimetric accuracies. Since more accurate models will be 
available when the satellite gradiometer mission is launched, then for 
solutions using this algorithm, the resulting accuracies will be 
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significantly better than those reported here. As gravity models 

improve, the degree and order of their expansions will increase. This 
can also improve the accuracies of the resulting gravimetric quantities 
computed by this algorithm, if the uncertanties of the coefficients 
remain optimally small. 

In summary, the primary investigation, in the form of a local 
collocation error analysis, has been described and complete results for a 
200 km altitude gradiometer mission were presented with and without the 
influence of the DEC 81 coefficient uncertainties. Th<» choice of degree 
variance model used in the study was made in lieu of available software 
and that the parameters of Jekeli adequately represent the actual 
gravity field in terms of the implied gravimetric variances. The effects 
of the DEC 81 geopotential coefficient uncertainties upon the degree 
variances were mentioned. It was noted that the error degree variances 

i 

abruptly changed the power spectrum from its smooth character implied 
by the modeled degree variances and further, the minimum case of 
assuming errorless coefficients was described as a useful method to 
determine the sensitivity of the resulting gravimetric accuracies. The 
accuracy results were presented and their characteristics were 
discussed particularly with regard to data density, mission altitude, and 
reference coefficient uncertainty sensitivity. In the final chapter, 
comparisons to other results will by made along with concluding remarks. 


5. DISCUSSION AND CONCLUSIONS 


5.1 Comparison to other Investigations 

It would be useful to place the results of this study in context with 
other similar investigations. Unfortunately, this is not easily done since 
the investigation by Reed [19721 and Krynski & Schwarz [1977] solve for 
different quantities than those determined here and since the initial 
assumptions of those investigation differ considerably from those made 
in this work. It is tempting to compare the results of this work with 
those of a similar study involving satellite to satellite tracking by Hajela 
[1983]. For 30'x30' areas, Hajela determined the accuracy of the 
predicted mean gravity anomaly to be 18.1 mgal. However, this result 
was for a mean gravity anomaly referred to the GEM 9 gravity field 
model. The results of the present report are referred to the 180 degree 
and order DEC 81 gravity model of Rapp. Hence, the accuracies 
reported in this work are expected to be better since the siodeled 
degree variances between degrees 20 and 180, as used in Hajela’s study, 
are replaced by the error degree variances implied by the DEC 81 
coefficient variances which have considerably smaller spectral power (cf 
Figure 11a). 

The rapid error analysis procedure of Jekeli A Rapp [1980] has been 
used in this report to generate Figures la, lb and 7. The procedure 
utilizes the Tscherning A Rapp degree variance model and for the 
construction of the figures mentioned above, Tscherning A Rapp's 
degree variance parameters have been used. A portion of the results 
used in the construction of the figures is given in Table 5 along with 
results from other runs of the procedure using Jekeli's parameters 
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(indicated by Jek). The procedure aasumea that the entire global data 
aet of meaaurementa in uaed in the computation of gravity anomaly or 
height anomaly accuracies. The accur*. iea computed by the procedure 
are expected to be better than thoae computed by the local collocation 
algorithm developed here. Thia ia especially true for the height anomaly 
accuracy since the procedure includes, in its computation, the long 
wavelength spectrum which comprise the bulk of the total height 
anomaly signal. The best accuracies from the local collocation algorithm 
tested in this study are (for 30'x30' means at gradiometer altitude of 
200 km, 10“ 3 E instrument accuracy): 9.1 mgal and 43.3 cm (Figures 13f 

and 13g) for gravity anomalies and height anomalies respectively, where 
the DEC 81 coefficient errors are included. The perfect gravity field 
case yielded best values of 8.9 mgal and 13.6 cm. From Table 5, the 
corresponding accuracies computed by the rapid analysis procedure are 
5.8 mgal and 9.4 cm. Therefore the differences of 3.1 mgal and 4.2 cm 
between the local collocation algorithm (with perfect potential 
coefficients) and the rapid analysis procedure is due to neglecting the 
region exterior to the grid in the local collocation algorithm as well as 
to the smoothing effects incured during the gridding of the data. The 
gravity anomaly accuracy can be improved by lowering the mission 
altitude, by increasing the instrument sensitivity, or by increasing the 
gridded data density. There are difficulties associated with all three 
ways to improve the accuracy. Lowering the altitude will increase the 
drag of the satellite thereby necessitating a disturbance compensation 
system similar to that proposed for GRM. Increasing instrument 
sensitivity may be possible with more research by instrument designers 
(some designers are optimistic that gradiometers may achieve 10~ s E 
sensitivity levels in the near future. The difficulties associated with 
atmospheric drag and gradiometer sensitivity are mission design 
considerations. The increasing of the gridded data density is more 
directly a data reduction problem. More specifically, when the gridded 
data is densified, a larger inversion process entails, thus increasing the 
computational demand. Two factors limit the gridded data density. One 
is physical, the other is economical. The physical limit is due to the 


73 


fact that tho data connot ba gridded more densely tUn tha a pacing of 
tha raw gradiant maaauramanta made in tha aatallita orbit. To do ao 
would introduce large errors in tha high frequency spectrum of tha 
gradianta. Tha economic limit ia one acceptable CPU-time required for 
the inversion. This of course, is a more subjective limit which depends 
on many factors. Nevertheless, it should be kept in mind during the 
selection of satellite grad iome ter data reduction algorithm. 


Table 5. 30'x30' mean accuracies implied by Jekeli k Rapp's rapid 

analysis procedure. 1 sec data integration period. 


Measurement 

Accuracy 

<B) 

Altitude 

(km) 

Param. 

DEC 81 
errors? 

TOTAL 
Qrav. Anon. 

(*«al> 

RMS ERRORS 
Height Anon, 
(cm) 

0 001 

160 

T/R 

N 

4.8 

6.3 

0.001 

160 

JBK 

N 

3.8 

5.2 

0.001 

200 

T/R 

N 

7.3 

11.4 

0.001 

200 

T/R 

Y 

7.5 

11.7 

0.001 

200 

JBK 

N 

5.8 

9.4 

0.0001 

160 

T/R 

N 

2.9 

3.1 

0.0001 

160 

JBK 

N 

2.2 

2.5 

0.0001 

200 

T/R 

N 

5.1 

6.7 

0.0001 

200 

JBK 

N 

3.9 

15.4 

*1.0 p/s 

160 

T/R 

N 

8.9 

15.6 

*1.0 p/s 

160 

JBK 

N 

7.1 

13.0 


*GHM parameters, 4s data integration, 300 km satellite separation 


5.2 Conclusions and Recommendations 

This report has described the data reduction of satellite born 
gradiomotry measurements by the local colloctMon method. Only the 

error analysis aspect has been investigated in this study since no test 
data yet exists. The algorithm was simplified to utilize only the 
anomalous vertical gravity gradients gridded at altitude. The anomalous 



gradient* are determined by removing reference gravity gradients 
computed from a high degree and order geopotential model. The 
uncertainties of the reference geopotential coefficients must be included 
in the covariance computations. Two cases were considered: one, where 

the uncertainties of Rapp's DEC 81 geopotential model [Fapp (1931)] were 
included in the covariance computations, and the other, where the 
gravity model to degree and order 180 was considered perfect. This 
was done to test the sensitivity of the resulting gravimetric accuracies 
to the uncertainties in the reference geopotential coefficients. It was 
shown that the height anomaly accuracy is indeed very sensitive to the 
influence of reference coefficient uncertainties which follows quite 
naturally since the height anomaly is a low frequency phenomenon and 
the coefficient uncertainties will contain these low frequencies. The 
gravity anomaly accuracy, on the other hand, did not improve as 
dramatically when the coefficient errors were removed. The results of 
this sensitivity study indicate that the accuracies of the gravimetric 
quantities computed by the local collocation algorithm significantly 
improve with better geopotential models. More accurate geopotential 
models should be available in the future when the satellite gradiometer 
mission is attempted. These models can be easily applied to the 
algorithm to provide more accurate gravimetric reductions. 

Further improvement may result with the inclusion of other 
gradient components (notably, those with at least one radial derivative). 
Krynski 4 Schwarz 11977] reported 15X to 20X accuracy improvement in 
the geoid undulation when the T , 3 and T a3 components were included. 
This could not be confirmed in this study due to software related 
difficulties. A drawback caused by inclusion of the T » 3 and T a3 
components results in the increased dimension of C 44 by a factor of 
three, further exacerbating the CPU time required for the matrix 
inversion. 

The difficulties with large matrix inversion has been mentioned 
several times in this report. Increasing the amount of data, either by 
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increasing data densities or by including other gradient components, will 
moat certainly increase the CPU time aa well as possibly degrading the 
stability of the inversion. Although these two aspects were not studied 
specifically, it is possible to suggest techniques to overcome these 
difficulties. The technique of frequency domain collocation provides a 
means for handling large amounts of data by reducing the computational 
demands incurred by the inversion process. The early study by Tait 
[1983] examined this possibility for airborne gravity gradiometry with 
rather promising results. This concept could be extended to satellite 
altitudes to test the validity of the technique. To stabilize the 
inversion process a new technique described by Jekeli [1985] may be 
very useful. The technique is known as Virtual Optimal Estimation 
which provides a stable method to invert an ill conditioned (or "almost" 
ill conditioned) symmetric definite matrix through iterative techniques. 
With thiB method, error introduced by an unstable inversion process can 
be avoided thereby strengthening the collocation solution. 

The local collocation algorithm was devised as an alternative method 
to compute surface mean anomalies based on data acquired by a satellite 
gradiometer mission. The accuracy results and ensuing discussion 
illustrates the behavior of the system with regard to data coverage, 
data density, altitude, and reference coefficient uncertainties. All of 
these factors play a crucial role in the accuracy of the computed mean 
gravity or height anomalies. To bring the mean gravity anomaly 
accuracies below the 5 mgal level using this algorithm, it is 
recommended that the mission altitude be kept as low as possible (e.g. 
160 km), the data be gridded as densely aB possible (e.g. <15‘), the 
gradiometer designed to be as sensitive as possible (e.g. <10~* E), and 
that the gravity model to be used for the computations of the reference 
gradients have the ssiallest coefficient uncertainties as possible. More 
research is needed to further refine the algorithm and to fully test its 
effectiveness by means of a simulation. 
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APPENDIX A 

Spherical Harmonic Expansions of the Reference Gradients 

In section (III. 2), it was shown that the "observed" gravity 

gradients in the local level coordinate system need to be centered in 
order to apply the least-squares collocation technique. The method to 
center the observations involves computing reference gravity gradients 
at points on the same geographical grid for which the observed 

gradients are known. The centered observation is then the difference 
of the observed gradient minus the reference gradient. In this 

appendix are found the expressions for the reference gradients given in 
terms of the reference geopotential model which are needed in equation 
(III.?). Reed [1973] originally derived the spherical harmonic 
expressions that ore given here for ease of reference. 

Using the local level coordinate system defined in Figure 3, the 

expressions are: 


N n 

Un(», X, r) = pj [-1 (f) n H (CnaCosmX + S^sinmX) P^i(sin»)] 

n=a m=o 

(A. 1) 

N n 

Uaa(*. X, r) = pr [“1 + YU (r)"!^ ( c na cos “X + S^ainmk) p££(sin*)] 

n=a a=o 

(A. 2) 
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N n 


U„(A, X, r)=SN[2 + r (friZ (C»co*jiX + S^inaX) P^CsinA)] 

n-a ■=« 

(A. 3) 

N n 


Uu(A» X, r) = pf ) (“) } (C^icoaX - S^sinaX) PnaCsinA)] 

1 2 «- ft 


n=a ■=© 


(A.4) 


w ii 

U. S (A, x, r) = pj } lf) n IZ ( c naCO««X - S,^sinaX) P^(sinA)] 


n=a a=o 


(A. 5) 


N n 


U a ,(A, X, r) = pf JZ (?)” 5ZZ ( c naCO**X + S^inaX) P?£(sinA)] 


n=a B=o 


(A. 6) 


where the N denotes the maximum degree of the reference geopotential 
model. The observed gradients can be expressed by allowing N to 
approach infinity and by substituting the unnormalized (conventional) 
geopotential coefficients C n m» s nm - or the conventional reference 
geopotential coefficients C^ m , S^ m . The superscripts attached to the 
conventional associated Legendre functions denote the differentiations. 
Reed [1973] alRo provides these differentiated functions in terms of 
non-differentiated Legendre functions: 


PA4(sin») = “ (n+» )]PnB(zinA) - tan* P„,«+i (sin») (A.7) 

Pg4(sii,#) = " (n+t) 2 ]Pn»(sinA) + tanA P nfB+1 (sinA) (A. 8) 


P*4(sinA) 


(n+l)(n+2)P nB (sinA) 


(A. 9) 
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PAi(»in#) 


■Os-l)sin» , a 

cos** P„*(sin*) + P„ (B+1 («!„♦) 


(A. 10) 


Pn«(*inA) = - P„ m (sin*) 


(A. 11) 


PSi(sinA) = a(n+a)tan* P nm (.!„♦) - (n+a) p^, 


(A. 12) 


The P nmisini) and Pn^sin*) terms can be found from familiar 
recursion formulas such as (Ilk [1983]); 


(n «)P M (t) = (an- i )tP n _, >a (t) - (n+a-i )P n -a tB (i) 


(A. 13) 


(i-t a )*p niB (t) = (aa-OtPn.^t) . ( n ^- I )(n-» ft )( l _t.)»p (t) 


(A. 14) 
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APPENDIX B 

Gradient Covariance Expressions 


In this appendix, the explicit covariance expressions are given for 
reference. The expressions are nothing more than applications of the 
covariance propagation law (Moritz, [1980]) 


Cij(P,Q) = L?LjK(P,Q) 


(B.l) 


from which all gravimetric covariances can be related to the disturbing 
potential covariance K(P,Q) by the functional L ( and Lj applied at the 
point P and Q. The covariance function of the disturbing potential is 
only a function of the relative location of the points P and Q. On the 
sphere, this can be written as (Moritz, [1980]) 


K(P,Q) = K(r P , r q , f Pfl ) = K 


(B.2) 


where r P and r 0 are the geocentric radial distances of the points P and 
Q and f PQ is the spherical distance of P and Q given in terms of their 
respective co-latitude and longitude as (0 - w/a - ♦), 


cosf P Q = cosdpcosdg + ain9 P sin0QCOs(X Q -X P ) 


(B.3) 


Thus the covariance of the disturbing potential is isotropic and 
stationary. 


Using the functionals given in equations (III. 15) to (III. 18), the 
covariances are determined and presented for reference below. Those 
functionals with a prime attached denote that the functional is to be 
applied «t the point Q. Unprimed quantities refer to P. 


I. Autocovariances of the observation covariance matrix, Ctt 


Cov[Tj » (P) »T, t (Q)] = [ 




r a cos a # *X a 


tan# * .1 
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* 4 K 
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* 3 K 
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cov[t 33 (p),t 33 (q)] = jpjps 


(B.6) 
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f » _ «’ , tan* * 1 

r * « a K . tan*' *K 1 

lr a cost «X*» r a *XJ 

lr' a cos*' «X#»' r' a «X'J 


* 4 K . tan»' «»K 


Hr^coikoir «♦' r a r' a coa* *X«a»X' 


. tan» * 3 K . tanitani' «»K 

rV^cost' •X#X'«»' r a r' a #X#X' 
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c.v[r„(P).T„(Q)] * jj] 
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i **K 
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II. Cross covariances of the observation covariance aatrix, Ctt* 
Cov[T„<P>.T„«») = r , r .,‘ cot ,, 

_ tani * 3 K _ tani * a K i_ * 3 K . i « a K ln . 

r a r' a •♦•♦' a r 2 r' «»*r rr' a *r**' a rr' *r*r' (B.10) 

Co»[T, 1 (P).T„(0)l = _Ii|. + i -£S_ (B.ll) 


Cov[f 1 1 (P) |T, a (Q) ] = r i r * a co8 a^ c08 ^' * r a r**eM*t *X a Jx 


tan#' « 3 K 


tan# **K tan# tan#' « a K 


r a r a 'coa#' *#*X'»#' r a r' a «##X' + rr' a coB#' #r#X'*#' 


• 3 K 


, tan#' * a K 
rr' a «r«X' 


(B. 12) 


Cov[T, i (P) |T, 3 (Q) ] - r a rc08 a* c08 *' #A a *x*«r' ~ r a r' a coa a #coa#' «A a *A' 


« 3 K 


tan# 


* 3 K tan# * a K i 

* T i * ’ » 9 x * T» T7 " ' 


« 3 K 


r a r'coa#' *#*X'*#' r a r' a cos#' *#*A' rr'coij' *r*X'*r' 


i * a K 

rr' a cos#' *r*A' 


(B. 13) 


Co»[T,i(P),I a ,(0)l = WIT - 7$T 


**K 


tan# * 3 K tan# * a K _j « 3 K 

i' a 3 TTTTT + 


1_ * a K 

r a r' *#*#'*r' T r a r' a *#<#' T rr' *r*#'*r' rr' a *r*#' 


(B. 14) 


Cov(T„(P),T„(0)) = JJTTpr * j 


(B. 15) 


Cov[T aa (P),T la (Q)] = r a r 'a cos ^' a*a # J' a #' + J^a # $aJ[ X ' 




* 3 K , tan#' * a K 

' m±' f * ! 


rr' a cos# <r*X'*#' rr' a 


*r«A 


(B.16) 
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Cov[T aa (P),T l3 (Q)} 


i » 

r a r'coa#' r a r' a coa#' *# a *X' 


. _» *fJC » * a K 

rr'coa#' *r«X'*r' rr' a coa# «r«X' 


(B. 17) 


Cov[I aa (P),T aa (Q)J 


» « 4 K i * 3 K . » * 3 K 

r a r' •r' “ r a r' a rr' *r' 



• »K 

•nr 


(B. 18) 


Cov[T„(P),T ia (Q)] = 


» * 4 K 

r' a coa#' *r a *X'*»' 


ton#' * 3 K 
r' a «r a *X' 


(B. 19) 


Cov[T aa (P),T ia (Q)J 


i i * 3 K 

r'coa#' *r a *X'*r' r' a coa#' *r a *X' 
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Cov[T 33 (P),T a3 (Q)] 


-L * 4 K 
r' #r a *#'*r r 


i « 3 K 
r' a *r a «#' 
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Cov[T ia (P),T l3 (Q)] 
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Cov[T ia (P),T a3 (Q)] 
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tan# « 3 K _ tan# « a K 

r a r' *X*#'*r' r a r' a *X##' 


(B.23) 


Cov[T 13 (P),T a3 (Q)] 


i « 4 K 

rr'cos# *X*r*#' a r' 


i » 3 K 

rr' a coa# *X*r*#' 


i * S K * a K 

r a r'coa# a X*#'*r' r a r' a co8# *X«#' 


(B.24) 
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III. Cross covariance matrix, C a t 


a. Gravity anosmlies, Ag. 


Cov[Ag(P) ,T, i (Q)) 


i tan»' * a K _ i_ * a K 

r' a cos a 4' *r*X' a r' a »r*4' r' *r*r' 


* a*K 2tan»' *K a «K 

rr' a cos a 4' *X' a rr' a *♦' rr' *r' 


(B.25) 


i #sv 

Cov[Ag(P),T aa (Q)] = - -ps 


_i a a K 
r' *r*r' 


a « a K a *K 
rr' a " rr' «r' 


(B.26) 


Covf Ag(P) ,T 33 (Q) ] 


* a K _ a * a K 
*r»r ' i 2 ~ r »r' a 


(B.27) 


Cov[Ag(P) ,T t a (Q) ] 


i # 3 K f nn4 |>V a |2v 

r' a cos4 *r*X'*4' ' "r 7 " araX/ “ rr' a cos4' a'x'ag' 


a *K_ 
rr' a cos4' »X' 


(B.28) 


Cov[Ag(P),T ls (Q)] 


i a»K » * a K a a»R 

r'cos4' *r*X'ar' r' a cos4' arax rr'cos#' ax'ar' 


+ 


a 

rr' a cos4' 


aK 

ax 


(B.29) 


Cov[Ag(P) ,T as (Q) ] 


i_ * 3 K i * a K _ _J_ aag a_ ag_ 

r' arag'ar' " r' a a r ag' rr' a*'# r ' rr' a a*' 

(B.30) 


B. Height anomaly f (** N the geoid height) 


ORIGINAL F • \ 

OF POOR QUALITY «6 

cov[f (p),Tn(o)J = ^ (prspr IT* - 3r * r 2H <»- 3I > 

Cov[f(P),T„(Q)] = “ [p? 7T? ♦ p- Jp-] (B.32) 

Cov[f(P),T„(0)] = ^ (B.33) 

Oo,(f(P).I„(0)] = ^ 7 $* + jH <»•») 

Cov[((P) ,T, j(O) ] = ^ [p^ - pr^ $•] (B.35) 

Cov[f(P),T„(Q)] = ^ [jr - ^7 jf-J (B.36) 


IV. The auto covariances of the predicted signals, C Ba . 

Fro* Moritz [1980], page 108, 

Ck>v[Ag(P) ,Ag(Q) ] = t ♦ p^ + ;?p' + 7 pK (B.37) 

and for the height, anosuly, 

Cov[f(P),f(0)J = “gK. 


(B.38) 


Appendix ", Listing of program 6IFRAD 

c 

c 

C PROGRAM GlFRAOl A PROGRAM WHICH COMPOTES THE EAPECTEO ERROR 
C IN A MEAN GRAVITY OR HEIGHT ANOMALY BASED 

C UPON RADIAL/RADIAl GRAD 10ME TERV MEASUREMENTS 

C EVER A GRID CP FIXED SUE. 

C 

C WRITTEN BY* JOHN N. ROBBINS 
C OATE* FEBRUARY 1. 1985 

C VERSt FEBRUARY 17* 1985 

C RUN I FEBRUARY 2 0, 1985 

C 

C DEPARTMENT OF GECDET1C SCIENCE AMD SURVEYING 
C OHIO STATE UNIVERSITY. 

C 

C 

IMPLICIT REAL *8(A-H*0-2I* LOGICAL (L) 

COMMON /1NTERP/CCV(3,519) 

DIMENSION CXX( 225*225) *DIST(225 V 225)*CSX(2*225> 

DIMENSION DSX (225)*SCR (225 )*CSX1 (1 *225)*CSX2(1 *225) 

DIMENSION CC1 (225*1 )*CC2(225*1) .ANSI l).ANT(i) 

OIMENSICN CXXI (225*225 I 

C DATA CSS1 *0552/0.11413538300. 224. 738999600/ 

DATA LSS1*CSS2/1 .1613600*3.021260402/ 

DATA CLAT/40 .00/ 

C 

C REAO IN THE COVARIANCES. 

C 

C FIRST* THE CXX ( CSX REQU1RE0 COVARIANCES ARE REAO. 

C 

REA0I5 *«)|C0V(1*II*C0V(2*1)*C0VI3*I)* 1-1*5191 
C 

c 

c ENO OF CXX C CSX RECUIREO COVARIANCES 
C 
C 
C 

C WRITE THE HEAOINC FOR THE GUTPUT TABLE. 

C 

C 

WRITE 16 *106 I 

106 FORMAT (2X* •GRID* *6X * *CR ID • *4X • 'MEAN GRAVITY ' *2X* *HEAN HGT *,3X, 
**1MV£RS.*»6X »*CSX*CXXI*CXS VALUES**/ 

♦2X • ‘ WIDTH* *3X * • SPAC INC • *5X *• ANOMALY* * 6X* • ANOMALY* *4X *• STA8IL. * •/ 
•2X* • fDEGI* *4X . *1 DEG ) • ,6X * * IMGAL) • »7X * * IME TERS9 • *1 2X**f METERS662) 
•8X* • IMGAL662I * *//l 
C 

C READ IN THE OVERALL SUE OF THE GRID ANO WHETHER OR NOT IT 

C IS THE LAST GRID SUE TO BE CONSIDERED (LOGICAL VARIABLE > LI. 

C 

3 REA0I5 *99)0 *OELO *L 

99 FORMAT (2F]0.7*L2I 

C WRITE (6 *98)0 *OELO*L 

Cfi FORMAT (//IX *2F 10.7*3X*L2) 

DDEC-0 

ODEL-OELD 

NP-I0NINT(0/0EL0H1 

NPS«NP*NP 

NPXl-(NPS-»l)/2 


«s- • 


NP12>NPI1«1 

INI UAL Ilf AU ARRAYS 10 UNO 10 AVOIO CARRY-OVER FROM 
PREVIOUS COMPILATIONS. 


ANSI!) 
AN? Ill 
DO 12 
00 12 
0 f ST f I 

cxxiii 

CXXII, 
00 1A 
00 1 A 

csxll. 

00 13 

CSX1I1 

CSX2I1 

CC1 ( I • 

CC2II. 

OSX(l) 


> 0.00 
• 0.00 
1>1 .225 
J-l .22 5 
•Jl-0.00 
• J) >0 .00 
J)-0 .00 
l>lr2 
J-I .225 
Jl-0.00 
1-1.225 

• n-o.oo 

• I) >0.00 
1 )> 0.00 
i)>o.oo 
-o.oo 


COMPUli INI MATRIX OF SPHERICAL DISTANCES FOR THE GR10 
SPECIFIED FOR THE OVERALL GRID SHE (0) AT THE LATITUDE 
CCLAT). THE SUBROUTINE ONLY NORMS WITH RECTANGULAR CRIOS. 
RAOIALLV SYMMETRIC CRIOS FUST BE DEALT WITH DIFFERENTLY. 

CALL S01SICLAT .O.NP.NPS.O I ST I 

SET-UP THE COVARIANCE MATRIX. CXX. 

00 11 1-1 .NPS 
00 11 J-I. NPS 

CXXII. J)*C0VIN1(01STf 1.JI.1) 

00 23 1-1 .NPS 
00 23 J-l.NPS 
CXXIJ.I)-CXX(I.J) 

A001TI0N OF INSTRUMENT NOISE TO DIAGONAL ELEMENTS. 

THE VALUE 1.0-6 IS THE SOUARE OF THE INSTRUMENT SENSITIVITY 
OF 1.0-3 EOTVOS UNITS. 

ONF. MUST REPLACE THIS VALUE FOR OTHER SENSITIVITIES. 

00 2A 1-1 .NPS 

cxxu.n-cxxf i.im.o-6 

SET-UP OF CSX RELATEO DISTANCES. 


00 21 I-1.NPX1 
DSXID-OlSTII .NPX1) 
00 22 I-NPX2.NPS 
OSXID-OISTINPXI .1) 


SET-UP THE CSX-MATRIX. 


00 31 J-l.NPS 

CSXf l.Jl-COVINTf OSX I J) .2) 


nnnftnonftnnrjnnnftooonnnnftnnftn non nr^onnooooon 
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31 CSX (2»«J) c C0V IN 1( OSX ( J) ,3) 

00 33 J-l.NPS 
CSXltl.J)-CSX(l.J) 

33 CSX2C1 . J)*CSX(2»J) 

HAKE CONFUTATIONS FOR THE EXPECTEO ERRORS RESULTING PROP 
THE GIVEN PARAMETERS. 


THE SUBROUTINES LIN V1F .VNLLFP, VMULFF ARE PR0V1DE0 BY 
IM5L SUBROUTINE LIBRARY. 

LINV1F * MATRIX INVERSION 

VHULFP, VMULFF * MATRIX ALGEBRAIC OPERATIONS. 

CALL LlNVlFCCXX,NPS.225,CX*I,f,StR,IER) 

CALL VMULFPICXXI ,CSX1 .NPS .HPS.l .225, 1 ,CCI ,223, 1ER1 ) 

CALL VMULF"(CXX1,CSX2,NPS,NPS,1,225,1,CC2,225,1ER2> 

CALL VMULFF (CSX1 ,CC1 ,1 ,NPS,1 ,1 ,225,ANS,1, IER3I 
CALL VMULFF CCSX2 ,CC2,1 ,NPS,1 ,1 ,223,AN1 ,1 ,IER*> 

FHA-CSSl-ANS(l) 

ECA-CSS2-ANTC I) 

FGA*OSGRT(EGAI 
FHA ■OSGRTfEHA ) 

WRITE RESULTS. 

URITE(6,10B)OOEG,OOEL,FGA,FHA,IER.ANSI1>,ANT(1) 

08 FORMAT IF6.5.3X ,F7.5,3X ,F10.5,3X,F3.5,4X,I4,3X,014 .7, 3X ,013.71 

IFI.NOT.LIGO TC 3 
STOP 
END 

SUBROUTINE SOI StCLAT ,0 ,NP ,NPS,0 1ST ) 


«*****«*«»*«*»«****««**#«*•«******+*******«»***•**#*** 

SUBROUTINE NOTES: VARIABLE LIST. 

1. INPUT VARIABLES IN THE CALL STATEMENT 9 

CLAT X THE LATITUOE CF THE CENTRAL POINT OF THE GRID. 

0 t THE OVERALL SI2E OF THE GRID. 

NP : THE NUMFER OF INCREMENTS ALONG A GRIO SIOE. 

NPS > THE TOTAL NUMBER OF POINTS WITHIN THE GRIO. 

2. OUTPUT VARIABLE RETURNEO TO MAIN PROGRAM! 

OIST : THE ARRAY OF SPHERICAL DISTANCES IN RADIANS,. 

IT IS UPPER” T$ ! ANGULAR FORM. 

3. VARIABLES USED WITHIN THE SUBROUTINE: 

NPM1 ! HP - 1 

ONP X 1./NPH1 - l./INP-ll, THIS GIVES THE RELATIVE 
GRIO SPACING VALUE. 


C 





NP2 * NPANPHld • H91*2 - NP ♦ 1 

NPMA > NPS - 1. TOTAL NO. OF POINIS MINUS 1. 

'***»**#*********#***#*« 4 **********«*********** 


IMPLICIT REAL «8U-H.0-2) 

DIMENSION DIS1 (2 25 $229 ) *PLAT|225).PLON(225) 

SD<N.X.V»2)-DACOS(DSIN<H)*OSINm«OCOS<M>*DCOS<VMDC03(2-X)) 

PI-4.DCHAT AN (1.001 

0«D*PI/1BO.OO 

BLAT-CLA1«P I/I80 .00 

NPH1-NP-1 

0NP-1.D0/NPM1 

NP2-NP4NFH141 

NPHA-NPS-1 

ESTADLISN LATITUDES AND LCNG1TU0ES FOR EACH OF THE POINIS. 

00 5 J-l *NP2 *NP 
JA-J-i 

0 JNP-FLCAT I JA 1 /FLOAT INP 1 
00 5 I-l.NP 

PLAT tI4JAI-BLAT40/2.00-OJNP*DNP4D 
00 6 I-l.NP 
I A- I-l 

00 6 J-1.NP2.NP 
JA-J-1 

PLONf 1 4JAI-FL0A.7 ( I A )*0NP*0 

COMPUTE THE SPHERICAL OISTANCES. 

00 11 I-l.NPMA 
II 

00 11 J-II.NPS 

DISTII.Jl-SDfPLATfl I.PLONI ll.PLAUJ) .PLONf JI) 

00 12 1*1. NPS 
D1ST(I .1 1-0.00 
RETURN 
END 

FUNCTION C0VIN1ID.I0) 

FUNCTION COVINT* INTERPOL A 1ES THE COVARIANCES REQUIREO IN 
THE SOLUTION. THE COVARIANCE TAPLES ARE PREVIOUSLY 
COMPUTED FROM THE *TORSICN BALANCE VERSION* OF COVAX 
BV C.C. TSCHERN1NC AND TABULATED FOR READING INTO THIS 
PROGRAM. THE FUNCTION HA r AS INPUT 1HE SPHERICAL 01STANCE 
•0* ANO THE COVARIANCE WE. MB'. IN THIS VERSION MO* 
HAS THE FOLLCHING MEANINGS* 

10-1 I RA01AL/RA0IAL CRA01ENT AUTC -COVARIANCES 
ID-2 t GRADIENT/HE IGHT ANOMALY CROSS-COVARIANCES 
ID-3 * GRADIENT /GRAVITY ANOMALY CROSS -COVARIANCES 

THE OUTPUT IS 'COVINT* MH1CH RETURNS TO THE MAIN PROGRAM. 

IMPLICIT REAL *8(A-H.0-2I 
COMMON /INTERP/C0VI3.S19) 
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\TT \ — •’ 


00 * 0.00 

0H-1.0-6 

IF (lO.EG.O)GO TO IS 
0A-0A6S10) 

IF (OA.LT.OM)CC TO 16 
IF (10.LT.0) CO TO 20 
NK-0 
GO TO 21 

20 HK-1 
10 - 10 * 1 - 1 ) 

21 OMIR-2 .9066208 70-4 
DR-O/OMIR 
X*D|NT(OR) 

F ft AC -OR -X 
IX-IOINT(OR) 

C 

C IN THE NEXT STATEMENT • THE VALUE *519* CORRESPONDS TO THE 

C MAXIMUM SPHERICAL DISTANCE (IN THIS CASE. 8.5 0E6REES) . 

C 

IF (IX.CT.S19) CO TO 17 
IX1*IX«) 

IX2*1X«2 

COVINT-CCV (ID. IXIMFRACMCCVC 10, IX2)-C0VC 10,1X1)) 

IF (MK.EG.l)COVINT— COVINT 
CO TO 16 

15 C0V1NT-0.00 

16 RETURN 

18 COVINT-COV(IO.l) 

RETURN 

17 HRlTE(6t96) ‘ 

98 FORMAT (///IX t 'THE SPHERICAL OISTANCE IS TOO LARGE*) f 

STOP 

ENO 


(* 
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